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An  understanding  and  quantitative  knowledge  of  the  flew  through  a  cascade 
of  airfoils  is  an  inport ant  component  in  the  gas  turbine  engine  design  process. 
Blade  passage  regions  in  both  turbine  and  compressor  stages  are  critical  areas 
in  which  improper  design  cay  lead  to  unacceptable  engine  performance.  There¬ 
fore  a  procedure  capable  of  predicting  the  flow  field  in  blade  passages  is  an 
important  consideration  for  the  design  engineer.  One  step  in  achieving  this 
objective  is  the  development  of  a  tve-c ir.ensional  cascade  flow  field  analysis. 

At  the  present  time  there  are  several  possible  approaches  to  the  cascade 
flow  analysis  problem.  These  include  (i)  inviscid  analyses,  (ii)  inviscid 
analyses  with  viscous  corrections  and  (iii)  fully  viscous  analyses.  Although 
inviscid  analyses  (e.g.  Refs.  1-6)  can  give  accurate  predictions  of  the  blade 
pressure  distribution  in  some  cases,  they  are  limited  by  their  neglect  of 
viscous  phenomena  which  for  some  applications  nay  be  important.  The  limita¬ 
tions  result  from  (i)  the  necessity  of  assuming  airfoil  circulation  to 
obtain  a  unique  flow  field  prediction,  (ii)  neglect  of  viscous  displacement 
effects,  and  (iii)  the  inability  to  predict  losses  and  heat  transfer.  Consider¬ 
ing  the  first  of  these  aspects,  specification  of  the  airfoil  circulation  may  be 
both  reasonable  and  s traigh t f orvard  for  isolated  sharp  trailing  edge  airfoils 
at  high  Reynolds  number  operating  below  stall;  in  these  cases  the  Kut ta-Joukow’ski 
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condition,  which  requires  the  trailing  edge  point  to  be  a  stagnation  point,  has 
proved  to  be  an  effective  tool  in  setting  the  airfoil  circulation.  However, 
proper  specif icat ion  of  the  circulation  is  not  obvious  for  rounded  trailing 
edge  airfoils,  such  as  those  typically  found  in  turbine  cascades.  Specifica¬ 
tion  also  remains  a  problem  in  time-dependent  flows  or  in  flows  about  blades 
having  trailing  edge  separation.  In  all  these  cases,  specification  of  the 
airfoil  circulation  may  require  use  of  an  experimental  data  base. 

Inviscid  flow  analyses  also  are  limited  by  their  neglect  of  viscous  dis¬ 
placement  effects.  In  cases  where  the  boundary  layer  developing  on  the  blade 
surface  remains  thin,  flow  field  predictions  neglecting  boundary  layer  displace¬ 
ment  effects  may  be  quite  accurate.  However,  if  the  boundary  layer  does  not 
remain  thin,  the  actual  pressure  distribution  generated  in  the  flow  may  be 
significantly  affected  by  the  boundary  layer  development.  In  these  cases, 
the  actual  pressure  distribution  is  often  well  represented  by  the  inviscid 
pressure  distribution  about  an  equivalent  body  consisting  of  the  physical  body 
plus  the  viscous  displacement  surface.  The  viscous  displacement  effects  will 
significantly  alter  the  generated  pressure  field  if  the  viscous  displacement 
thickness  becomes  a  significant  fraction  of  the  blade  thickness,  as  is  the  case 
when  boundary  layer  separation  occurs.  Even  when  viscous  displacement  effects 
are  small,  they  may  still  be  important  and  a  particular  example  of  such  a 
situation  occurs  in  transonic  flow.  In  transonic  flow  the  flow  field  is 
extremely  sensitive  to  the  effective  airfoil  shape,  and  an  inviscid  determin¬ 
ation  of  the  shock  location  may  be  unreliable  unless  the  boundary  layer  thick¬ 
ness  is  very  small.  When  the  physical  boundary  layer  viscous  displacement 
effects  change  the  effective  airfoil  geometry,  an  inviscid  analysis  prediction 
of  shock  location  and  strength  may  be  significantly  in  error  leading  to  a 
poor  prediction  of  the  entire  flow  field.  A  final  problem  associated  with 
inviscid  analyses  is  their  inability  to  predict  aerodynamic  losses  and  heat 
transfer  rates.  If  these  quantities  are  desired,  the  inviscid  analysis  must 
be  combined  with  some  viscous  correction  procedure. 

The  need  for  including  viscous  phenomena  in  a  cascade  analysis  has  led  to 
several  approaches.  One  approach  modifies  the  inviscid  solutions  for  viscous 
effects  via  empirical  data  correlations  or  solves  the  boundary  layer  equations 
under  the  calculated  inviscid  pressure  distribution.  In  regard  to  the  first 
of  these,  although  methods  based  upon  data  correlations  can  play  a  useful  role 
in  the  design  process,  they  are  limited  to  flow  conditions  within  the  range  of 
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correlating  data.  The  second  approach  is  considerably  ncre  gcr.eral  and  a 
review  of  boundary  layer  techniques  for  turbomachinery  applications  has  been 
presented  by  McDonald  (Ref.  7).  When  the  boundary  layer  retrains  thin  through¬ 
out  the  flow,  as  it  does  in  many  practical  applications,  the  combined  inviscid 
flow-boundary  layer  analysis  may  give  an  adequate  description  of  the  flow  even 
if  viscous  displacement  effects  on  the  pressure  distribution  are  neglected. 
However,  as  previously  discussed,  in  many  practical  applications  severe 
pressure  gradients,  shock  waves  or  transonic  effects  can  cause  significant 
displacement  effects  which  in  turn  result  in  the  actual  pressure  distribution 
being  considerably  different  from  that  obtained  via  a  purely  inviscid  analysis. 
In  these  situations  an  analysis  which  solves  an  inviscid  set  of  equations  to 
obtain  a  blade  pressure  distribution  and  then  solves  a  set  of  blade  boundary 
layer  equations  under  the  influence  of  this  pressure  distribution  to  account 
for  viscous  effects  is  inadequate,  and  a  solution  recognizing  the  mutual 
dependence  of  the  pressure  distribution  and  blade  viscous  effects  is  required. 
Two  major  approaches  have  been  developed  for  such  problems;  these  are  the 
inviscid  flow-boundary  layer  strong  interaction  analysis,  and  the  fully 
viscous  passage  analysis.  Rather  than  neglect  the  effect  of  viscous  displace¬ 
ment  upon  the  pressure  field,  boundary-layer  strong-interaction  analyses  solve 
for  both  the  nominally  inviscid  outer  flow  and  the  wall  viscous  layer  in  such 
a  manner  that  viscous  displacement  effects  are  allowed  to  influence  the  gener¬ 
ated  pressure  distribution.  Typical  strong  interaction  procedures  can  be 
found  in  Refs.  8-10. 

A  strong  interaction  analysis  may  take  the  form  of  either  a  forward¬ 
marching  procedure  or  a  global  iteration.  For  regions  where  the  inviscid  flow 
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is  supersonic  and  thus  described  by  hyperbolic  equations,  a  solution  can  be 
f orward-marcned  in  space  with  the  inviscid  and  viscous  regions  coupled  on  a 
station-bv-station  basis.  The  chief  difficulty  in  this  process  is  that  math¬ 
ematically  stiff  equations  must  be  solved.  Common  problems  with  stiff  equations 
show  up  in  the  form  of  numerical  solutions  which  can  quickly  branch  off  the 
desired  solution  thus  producing  a  physically  unrealistic  result.  In  regions 
where  the  inviscid  flow  is  subsonic  and  thus  described  by  elliptic  equations, 
a  forward-marching  procedure  without  iteration  is  not  physically  realistic  as 
a  result  of  the  upstream  pressure  propagation,  and  consequently  a  sequence  of 
inviscid  and  boundary  layer  solutions  must  be  performed  in  a  manner  where  each 
stage  corrects  the  former  one  through  a  global  iteration.  Additional  problems 
occur  in  transonic  cascades  where  mixed  regions  of  flow  are  present  and  where 
small  displacement  effects  may  have  considerable  influence  on  shock  location 
and  overall  pressure  distribution.  Further,  if  an  interaction  procedure  is  to 
be  used,  the  viscous  layer  can  be  solved  iteratively  by  either  a  forward¬ 
marching  boundary  layer  calculation  procedure  or  as  the  asymptotic  condition 
of  a  time  or  pseudo-time  dependent  integration.  In  the  case  of  steady  state 
forward-marching  boundary  layer  procedures,  problems  are  encountered  when  the 
boundary  layer  is  subjected  to  a  strong  enough  adverse  pressure  gradient  to 
cause  separation.  Although  a  boundary  layer  procedure  can  be  forward  marched 
through  separation  by  suppressing  streamwise  convection  terms  in  the  separated 
flow  region  (the  FLARE  approximation)  (Ref.  11),  the  resulting  solution  is 
based  upon  an  approximation  made  in  the  separated  flow  region  which  becomes 
progressively  more  inaccurate  as  the  backflow  velocities  increase.  Therefore, 
calculated  details  of  the  flow  in  this  region  cannot  be  expected  to  be  accurate 
with  significant  backflow  using  FLARE.  A  global,  but  time  consuming,  iteration 
is  necessary  to  replace  the  FLARE  approximat ion ,  and  at  this  point  the  effi¬ 
ciency  of  the  forward-marching  scheme  must  be  carefully  evaluated  to  ensure  a 
net  gain  exists  relative  to  solving  the  full  Navier-Stokes  equations  for  the 
viscous  layer.  Time  integration  of  the  interaction  between  the  boundary  layer 
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and  the  external  flow  can  be  structured  to  avoid  the  use  of  either  F! A*- F  r 
the  global  iteration  to  account  specifically  for  the  backflow  velocity  f  F-_ :  . 
and  in  some  cases  the  time  marched  iteration  approach  may  show*  a  significant 
gain  relative  to  solving  the  Navier-Stokes  equations  for  the  viscous  flow 
region.  However,  the  interaction  approach  remains  limited,  even  with  tine 
integration,  to  flows  where  the  division  of  the  flow  into  zones  which  can  be 
interacted  is  reasonable. 

The  second  approach  used  in  obtaining  solutions  for  flow's  in  which  viscous 
phenomena  significantly  affect  pressure  distribution  is  the  fully  viscous 
analysis.  The  fully  viscous  analysis  solves  the  entire  flow  field  via  a  single 
set  of  equations,  thus  avoiding  any  need  to  divide  the  flow  into  viscous  and 
inviscid  regions.  Two  general  approaches  of  this  type  for  cascades  have  been 
used;  one  approach  by  Shararoth,  Gibeling  and  McDonald  (Ref.  13),  is  based 
upon  the  full  ensemble  averaged  Navier-Stokes  equations  and  a  second  approach 
based  upon  the  so-called  thin  shear  layer  equations  has  been  used  by  Steger, 
Pulliam  and  Chiraa  (Ref.  6).  The  thin  shear  layer  equations  are  an  approximate 
form  of  the  Navier-Stokes  equations  which  contain  all  pressure  and  convective 
terms  but  retain  only  those  viscous  terms  significant  in  thin  shear  layer  flow 
aligned  with  one  coordinate;  the  remaining  viscous  terms  are  omitted  from  the 
analysis.  The  fully  viscous  passage  analyses  overcome  many  of  the  problems 
associated  with  boundary  layer  strong  interaction  analyses  as  they  avoid  the 
need  to  divide  the  flow  into  viscous  and  inviscid  regions  and  eliminate  the 
iterative  matching  between  regions.  Furthermore,  the  full  Navier-Stokes 
analysis  contains  all  necessary  convective,  pressure  and  diffusion  terms  neces¬ 
sary  to  describe  large  separated  regions  correctly  and,  therefore,  flows  with 
significant  separation  or  shear  layers  not  aligned  with  one  of  the  coordinates 
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present  no  special  problem  when  the  full  set  of  equations  is  used.  Further¬ 
more,  since  fully  viscous  passage  analyses  are  based  upon  tine-dependent  forms 
of  the  governing  equations,  these  approaches  can  be  used  to  study  steady  state 
flows,  should  they  exist,  as  the  asymptote  of  time  integration  and,  in  addition, 
these  procedures  can  be  used  to  study  transient  effects. 

Demonstrat ions  of  thin  shear  layer  calculations  are  presented  in  Ref.  6, 
and  demonstrations  of  Navier-Stokes  calculations  are  presented  in  Ref.  13.  The 
results  of  Ref.  13  show  that  for  this  problem,  the  Linearized  Block  Implicit  (LBI) 
procedure  of  Briley  &  McDonald  (Ref.  14)  was  very  efficient  in  solving  the  full 
Navier-Stokes  equations,  and  converged  pressure  distributions  for  an  unstaggered 
cascade  of  NACA  0012  airfoils  were  obtained  in  low  Mach  number  laminar  flow 
within  sixty  time  steps.  Turbulent  flow  calculations  for  the  ensemble-averaged 
Navier-Stokes  equations  also  converged  within  sixty  time  steps,  with  the  cal¬ 
culation  being  initiated  from  a  converged  laminar  solution.  These  very  fast 
convergence  rates  demonstrate  the  practicality  of  applying  Navier-Stokes 
analyses  to  the  cascade  problem.  The  results  presented  in  Ref.  13  considered 
subsonic  laminar  and  turbulent  flow  (M^  ~  0.15)  in  a  very  simple  cascade  of 
unstaggered  NACA0012  airfoils.  The  present  report  extends  these  results  to 
transonic  flow  in  a  realistic  cascade  geometry.  As  a  prelude  to  the  transonic 
cascade  calculation,  it  was  deemed  appropriate  to  consider  a  simple  one- 
dimensional  transonic  flow  and  use  this  simple  transonic  test  case  to  assess 
the  role  of  type-dependent  differencing,  boundary  conditions  and  artificial 
dissipation  techniques.  The  present  report  discusses  this  assessment  as  well 
as  extension  of  the  constructive  coordinate  system  generation  process  to 
include  realistic  cascade  configurations.  Finally,  a  sample  high  Reynolds 
number  transonic  calculation  is  presented. 


14.  Briley,  W.R.  and  H.  McDonald:  Solution  of  the  Mil t idimens ional  Compres¬ 
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Physics,  Vol.  24,  No.  4,  August  1977,  p.  372. 


II  NUMERICAL  CONSIDERATIONS 
Tine-Marching  and  Iterative  cr  Relaxation  Pro.- ---cures 

To  date  both  tine-marching  and  iterative  or  relaxation  solutions  have  been 
used  to  obtain  steady  solutions  to  transonic  flow  problems.  In  the  tine  march¬ 
ing  approach  the  governing  unsteady  equations  are  solved  under  the  influence  of 
steady  state  boundary  conditions  and  a  steady  flow  field,  if  it  exists,  may  de¬ 
velop  as  the  solution  is  marched  in  tine  (e.g.  Refs.  13  -  13).  In  the  relaxa¬ 
tion  approach  the  steady  state  form  of  the  equations  is  considered  and  a 
converged  solution,  if  it  can  be  obtained,  is  achieved  by  first  assuming  values 
for  the  dependent  variables  throughout  the  flow  field  and  then  iteratively  cor¬ 
recting  or  "relaxing"  the  assumed  solution  in  a  series  of  iterations  until  the 
governing  steady  state  equations  are  satisfied  (e.g.  Refs.  16  -  19).  Although 
each  iterative  correction  or  relaxation  step  can  often  be  considered  in  some 
sense  as  a  pseudo  time  step  (e.g.  Ref.  16),  the  time  marching  and  iterative 
methods  are  not  necessarily  equivalent.  One  obvious  difference  concerns  cases 
where  no  steady  flow  exists.  In  these  cases  an  iterative  or  relaxation  proce¬ 
dure,  unless  it  can  be  formally  related  to  the  transient  evolution,  can  at  best 
only  give  some  mean  flow  upon  convergence  and  cannot  predict  the  experimentally 
observed  unsteady  flow  field.  However,  even  in  those  cases  where  a  steady  flow 
exists,  it  should  be  recognized  that  the  process  by  which  an  iterative  or  relaxa¬ 
tion  procedure  reaches  a  steady  solution  need  not  (and,  in  general,  does  not) 
conform  to  any  physical  process.  As  a  result,  physical  reasoning  based  upon 
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ritual  flow  processes  nay  not  be  a  valid  guide  in  developing  a  convergent  itera¬ 
tive  procedure.  The  converse  is  also  a  possibility,  that  strategics  and  argue- 
nents  developed  from  relaxation  approach  need  not  be  required  or  useful  in 
developing  the  transient  approach.  As  a  final  introductory  jonnent ,  it  is  useiu 
to  consider  two  subsets  of  time  marching  methods.  In  the  first  of  these,  no 
attempt  is  made  to  follow  transients  and  as  a  consequence  this  subset  nav  be 
thought  of  as  a  time  consistent  iterative  approach  to  solve  the  steady  flow 
equations.  In  the  second  subset,  again  no  attempt  is  made  to  follow  the  tran¬ 
sient  but  in  addition  the  time  step  may  be  varied  spatially  and  temporally  to 
improve  convergence.  This  is  equivalent  to  a  matrix  conditioning  (Ref.  20), 
but  with  the  special  constraint  that  the  conditioning  not  change  the  type  of 
the  conditioned  system  and  consequently  the  resulting  iteration  remains  time 
consistent.  For  reasons  that  will  be  subsequently  evident,  the  critical  dis¬ 
tinction  between  iterative  and  time  consistent  categories  lies  in  the  type  of 
inferred  equation  system  treated  in  the  evolution  process. 

An  example  demonstrating  these  points  can  be  found  in  the  two-dimensional 
velocity  potential  equation.  For  unsteady  flow7,  the  two-dimensional  velocity 
potential  equation  can  be  written  as 


r  tt 
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As  can  be  demonstrated  by  a  characteristic  analysis,  Eq .  (1)  is  a  second-order 
hyperbolic  equation  in  the  time  variable,  t.  If  the  transient  as  well  as  the 
steady  state  solution  of  velocity  potential  is  desired,  it  would  be  necessary 
consider  the  time-dependent  equation  as  represented  by  Eq .  (1).  However,  if 

only  the  steady  state  solution  is  of  interest,  an  iterative  or  relaxation  pro¬ 
cedure  can  be  applied  to  the  steady  form  of  Eq .  (2)  which  is  obtained  by  setting 
time  derivatives  to  zero 

2  2  2  2 
0  =  (a  ~  $  )  4>  +  (a  -  $  )  4>  (2) 

u  x  rxx  x  y  xy  ry  ryy  v  y 

When  Eq .  (2)  is  solved  by  a  relaxation  procedure,  values  are  assumed  for 
<j>(x,y)  at  every  grid  point  and  the  values  are  then  corrected  through  a  series 
of  iterations  until  a  steady  state  is  reached.  As  a  simple  example,  if  a  point 
relaxation  based  upon  simultaneous  displacement  were  used  to  solve  the  steady 
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z\:z -d  ir.ens  ional  potential  equation,  then  :*  would  be  related  ro  rv  an 
equation  of  the  form 
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where  n  is  the  iteration  index,  A,  B,  C  and  E  are  nonlinear  coefficients  and  D 
is  a  function  of  the  nonlinear  coefficients,  the  spatial  grid  size  and  the 
relaxation  parameter.  If  D  is  written  instead  as  1/4t  and  it  =  cn  ^-cn,  the 
difference  equation  is  a  consistent  approxima t ion  to  the  differential  equation 


$■  =A$  +  B$ 
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which,  unlike  the  transient  potential  equation,  is  a  parabolic  differential 
equation  of  the  first  order  in  the  marching  direction  t.  Equation  (4)  is  termed 
the  ’inferred  differential  equation’  associated  with  the  relaxation  procedure. 
Obviously,  this  equation  need  not  correspond  to  any  physical  transient  process 
during  its  iterative  development.  For  example,  using  a  change  of  independent 
variables,  Eq.  (4),  can  be  written  in  the  form 


*  =Ar4>,  ,+C’4>,  t+E’ 

t  xx  y  y 


Based  upon  experience  with  parabolic  equations  such  as  the  heat  conduction  or 
diffusion  equation,  numerical  difficulties  may  be  expected  (and  have  been 
experienced)  if  either  A’  or  C’  is  negative.  Further,  based  upon  Eq .  (2), 
negative  diffusion  is  clearly  a  possibility  depending  on  precisely  how  the 
iteration  is  set  up,  whenever  the  local  velocity  is  supersonic.  Any  such  nega¬ 
tive  diffusion  in  iteration  space,  if  it  were  to  occur,  need  have  no  physical 
relevance  to  the  problem  of  interest.  As  a  collorary,  numetical  strategies 
developed  to  counter  the  problem  of  negative  numerical  diffusion  in  interation 
space,  need  have  no  counterpart  in  the  transient  approach. 

Although  point  relaxation  with  simultaneous  displacement  serves  to 
demonstrate  some  differences  between  iterative  and  time-marching  procedures, 
point  relaxation  is  not  a  current  method  of  solving  transonic  potential  flow 
problems.  Much  more  efficient  methods  have  developed  and  are  being  used 
extensively  which  do  not  attempt  to  retain  time  consistency  and  are  based 
upon  successive  line  over-relaxation  (Refs.  16  and  17)  or  alternating 
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Jirevii^D  implicit  techniques  (Refs.  18  and  19).  A  line  re  I  a::n:  i  :-r. 
which  selves  the  tvo-c imensional  steady  velocity  potential  equation  u. -.ru  type- 
dependent  rotated  differences  has  been  examined  in  derail  by  Jameson  iRef.  17). 
The  resulting  inferred  differential  equation  is  typical  of  these  nore  efficient 
and  icjre  modern  time  inconsistent  iterative  schemes.  As  an  example  the  inferred 
differential  equation  has  [17]  the  principal  part 


(M  -IX  +  2ao 
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where  s  is  the  streamvise  direction,  n  is  the  direction  normal  to  a  streamline 
and  x  is  a  pseudo-time  variable  representing  the  relaxation  parameter.  The 
coefficients  a  and  6  depend  upon  the  relaxation  parameter.  Once  again  the  dif¬ 
ferential  equation  inferred  by  the  relaxation  process  does  not  represent  the 
C ine-dependent  velocity  potential  equation,  and  need  not  represent  any  physical 
process  as  it  evolves.  In  order  to  elucidate  Eq .  (5),  a  new  independent  vari¬ 

able,  T,  is  introduced 


~  _  US  ,  ^ 
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which  transforms  the  principal  part  to  the  form 
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The  resulting  equation,  like  the  physical  time-dependent  equation,  Eq .  (1), 
and  unlike  the  point  relaxation  equation,  Eq .  (4),  contains  second  derivatives 

in  the  three  dependent  variables  s,  n  and  T.  In  subsonic  flow  (M<1)  ,  Eq .  (7) 

is  hyperbolic  with  T  being  the  time-like  coordinate  and  is  of  a  similar  type 
as  the  transient  potential  equation.  In  supersonic  flow  (M>1)  the  inferred 
differential  equation  remains  hyperbolic,  however,  the  time-like  direction  is 
no  longer  the  independent  variable  T  but  rather  is  the  independent  variable  s, 
the  physical  streainwise  direction.  It  follows,  therefore,  that  strategies  re¬ 
quired  by  these  relaxation  schemes  to  preserve  stability  when  supersonic  regions 
are  encountered  are  not  required  a  priori  for  the  transient  approach  which  does 
not  change  type  upon  encountering  a  supersonic  region. 
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Although  iterative  and  relaxation  techniques  have  te_rs  very  :  u  r  ful 
in  obtaining  inexpensive  and  accurate  solutions  for  the  transonic  :::-,n tial 
flow  equation,  several  points  should  be  emphasized.  First,  the  inf-rr^d  dif¬ 
ferential  equation  (or  pseudo-time  differential  equation)  actually  being 
solved  is  not  time  consistent  in  supersonic  flow.  Secondly,  the  solution  to 
this  equation  only  represents  the  desired  physical  process  when  and  if  a  steady 
solution  to  the  equation  is  obtained.  Finally,  since  the  differential  equation 
being  solved  does  not  correspond  to  a  physical  process  during  its  transient 
phase,  a  physical  knowledge  of  the  actual  flow  process  need  not  be  relevant  in 
developing  a  convergent  solution  procedure.  Guidance  for  developing  stable 
algorithms,  boundary  conditions,  etc.  must  be  obtained  solely  from  the  mathe¬ 
matical  properties  of  the  inferred  differential  equation,  which  in  turn  must 
be  derived  and  analyzed.  Unlike  the  iterative  or  relaxation  methods,  the 
methods  which  solve  the  physical  unsteady  equations  model  a  physical  problem 
during  the  transient  flow  process.  Consequently ,  an  unders tanding  of  both  the 
flow  physics  and  the  mathematical  properties  of  the  equations  can  be  used  for 
guidance  in  choosing  algorithms,  boundary  conditions,  etc. 

The  foregoing  discussion  has  served  to  demonstrate  some  of  the  broad 
similarities  and  differences  between  iterative  and  tine-marching  or  time-consis¬ 
tent  procedures.  Both  procedures  have  been  used  to  solve  transonic  flow 
problems  (Refs.  13  -  19)  and  in  either  case  the  key  is  to  utilize  a  procedure 
which  allows  the  efficient  attainment  of  a  stable  steady  solution  consistent 
with  the  steady  differential  equation.  Considering  first  the  requirement  of 
efficiency,  the  procedure  should  converge  to  the  steady  solution  as  quickly  as 
possible  so  as  to  minimize  computer  run  time.  If  it  is  desired  to  follow  the 
physical  transients,  then  a  time-marching  solution  must  be  used  w’ith  the  time 
step  chosen  so  as  to  obtain  both  stability  and  transient  accuracy.  However, 
if  only  the  steady  solution  is  of  interest,  then  either  a  time-consistent 
solution  or  an  iterative  procedure  can  be  used.  However,  it  should  be  noted 
that  if  a  time-marching  procedure  is  utilized  solely  to  obtain  the  steady  solu- 
ion,  there  is  no  need  to  follow  the  transient  accurately.  The  only  requirement 
is  the  rapid  attainment  of  an  accurate  steady  solution,  and  the  time  step  may  be 
as  large  as  allowed  by  numerical  stability  considerat ions  and,  in  fact,  should 
be  chosen  to  optimize  the  convergence  rate.  The  experience  of  Shamroth, 

Gibeling  and  McDonald  with  subsonic  flows  (Ref.  13)  has  shown  cascade  flow 
fields  to  converge  to  steady  state  in  approximately  60  time  steps  using  a  split 
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Linearized  Block  Implicit  (LBI )  (Kef.  1-)  scheme  for  the  unsteady  Luvier- 
Stokes  equations.  As  demonstrated  subsequently  in  the  present  report  and 
although  no  attempt  has  yet  been  made  to  accelerate  the  convergence  rate  for 
transonic  cascade  calculations  for  this  flow,  convergence  is  still  obtained 
within  150  time  steps. 

In  the  past,  the  statement  has  frequent lv  been  made  that  iterative  or 
relaxation  schemes  converge  to  the  steady  state  more  quickly  than  do  tine 
marching  schemes.  On  the  basis  of  the  more  recent  studies  (Ref.  14),  including 
the  present,  it  is  clear  that  this  statment  applied  to  previous  stability 
restricted  explicit  t ime-marching  schemes,  such  as  Ref.  15.  This  statement 
is  not,  at  present,  applicable  to  some  of  the  current  generation  of  implicit 
time-marching  schemes.  Further,  it  is  also  clear  that  rapid  achievement  of 
the  steady  state  by  a  t ime-marching  method  results  in,  and  perhaps  even  requires 
that,  the  transients  not  be  resolved.  When  transients  are  not  resolved  the 
time-marching  approach  may  be  thought  of  as  time-consistent  iteration  and  as 
such  retains  the  advantage  that  only  the  boundary  conditions  need  distinguish 
between  subsonic  and  supersonic  flow.  If  the  transients  are  not  computed  ac¬ 
curately,  then  the  computed  steady  solution  is  subject  to  the  same  existence 
and  uniqueness  questions  that  the  relaxation  solutions  are.  The  primary’  advan¬ 
tage  of  the  time-marching  or  time-consistent  approach  is  the  ready  ability  to 
compute  the  transient  flow,  in  addition  to  the  steady  flow,  if  required.  This 
additional  feature  is  obtained  at  little  computational  disadvantage,  if  any, 
compared  to  the  current  generation  of  iterative  or  relaxation  approaches.  It 
is  also  observed  that  the  transient  approach  is  straight f orward  to  implement, 
and  does  not  require  special  distinction  in  the  interior  flow  between  subsonic 
and  supersonic  regions.  Boundary  conditions  are  another  matter,  and  in  accord¬ 
ance  with  the  flow  of  information  as  determined  by  the  local  characteristics, 
a  distinction  must  be  made  between  supersonic  and  subsonic  flow  on  the  boundaries. 

It  is  noted  that  conditioning  of  the  matrix  which  arises  from  an  implicit 
time-dependent  formulation  can  be  used  to  accelerate  convergence  (Ref.  20). 

When  this  conditioning  is  restricted  to  preraul tiplicat ion  of  the  matrix  eaua- 


20.  McDonald,  H.  and  W.  R.  Briley:  Computational  Fluid  Dynamic  Aspects  of 
Internal  Flows.  ALAA  Paper  79-1445. 
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zisv.  bv  3  positive  diagonal  matrix,  the  result  is  t*u  nival  ext  r 
spatiallv  varying  time  step  and  the  scheme  remains  tir.tr  e  L.-r.t  .  It  is 

readily  apparent  that  such  a  procedure  will  scale  the  slope  cf  th*.*  ch.arac:  er¬ 
istics  in  tine,  but  leaves  both  the  sign  and  the  eventual  steady  flow  unaltered. 
Such  a  conditioning  thus  leaves  the  type  of  the  governing  equations  unchanged, 
together  with  the  boundary  conditions.  Hence,  this  scaling  retains  the  attri- 
butes  of  the  transient  approach,  but  with  improved  convergence,  w’hile  retaining 
the  ability  to  perform  accurate  transient  calculations  by  simply  preeona i t ion- 
ing  with  the  unit  matrix. 

Spatial  Differences,  Boundary  Conditions  and  Artificial  Dissipation 

The  previous  favorable  experience  obtained  in  applying  the  LBI  procedure 
of  Ref.  14  to  subsonic  cascade  flow  fields  (Ref.  13)  in  conjunction  with  the 
considerations  of  the  previous  subsection  argue  strongly  for  applying  this 
same  technique  to  the  transonic  cascade  problem.  However,  prior  to  this 
application  several  items  must  be  considered.  These  include  spatial  differenc¬ 
ing,  boundary  conditions  and  artificial  dissipation.  Obviously,  these  items 
are  strongly  connected  and  cannot  be  considered  independently.  For  example, 
some  spatial  difference  schemes  contain  a  dissipative  truncation  error  which 
is  equivalent  to  an  implicit  artificial  dissipation.  Also  the  type  of  spatial 
differencing  chosen  will  affect  the  number  of  boundary  conditions  required  to 
close  the  algebraic  set  of  difference  equations. 

Previous  investigators  of  transonic  flow  problems  with  some  exceptions 
(e.g.  Ref.  15)  have  concentrated  primarily  upon  solutions  of  the  potential 
flow  equations.  Their  analyses  have  used  a  variety  of  combinations  for  spatial 
differencing,  boundary  conditions  and  artificial  dissipation.  The  present 
approach  focuses  upon  the  full  ensemble-averaged,  time-dependent  Navier-Stokes 
equations,  thus  including  both  rotational  and  viscous  effects.  The  analysis 
applies  a  split  LBI  algorithm  to  the  governing  equations,  expresses  spatial 
derivatives  via  central  differences  throughout  the  flow  and  selectively  adds 
an  explicit  artificial  dissipation  term  to  the  equations.  As  shall  be  demonstra¬ 
ted  subsequently,  when  this  procedure  is  used  with  careful  attention  paid  to 
boundary  conditions,  it  leads  to  rapid  convergence  and  sharp  representation  of 
the  transonic  shock.  Prior  to  consideration  of  the  Navier-S tokes  cascade  cal¬ 
culation,  the  problems  of  boundary  conditions,  spatial  differences  and  artifi¬ 
cial  dissipation  were  investigated  through  a  one-dimensional  transonic 
problem. 
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Model  Problem  for  Tran*; 


A  simple  nodel  problem  consisting  of  one-dimensional  flow  with  specified 
heat  sources  is  used  to  examine  techniques  for  treating  spatial  differencing, 
boundary  conditions  and  artificial  dissipation.  The  equations  solved  are  the 
continuity  equation 
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the  momentum  equation 
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and  the  energy  equation 
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wfhere  p  is  density,  u  is  velocity,  p  is  pressure,  is  specific  heat,  T  is 
total  temperature,  t  is  shear  stress,  q  is  a  heat  source,  t  is  time  and  x  is 
the  spatial  coordinate.  In  one-dimensional  flow,  acceleration  or  deceleration 
can  be  obtained  by  varying  the  heat  source-sink  distribution  q(x).  At  subsonic 
Mach  numbers,  a  heat  source  (q>0)  accelerates  the  flow  while  a  heat  (q<0)  sink 
decelerates  the  flow;  at  supersonic  Mach  numbers  the  effects  are  reversed  with 
heat  addition  leading  to  deceleration  and  heat  removal  leading  to  acceleration 
(cf.  Ref.  21).  By  an  appropriate  selection  of  Mach  number,  Reynolds  number  and 
heat  source-sink  distribution,  this  model  problem  permits  study  of  a  vide  range 
of  flow  conditions  including  subsonic,  supersonic  and  transonic  flow  w’ith  and 
without  discontinuities . 


21.  Shapiro,  A.:  The  Dynamics  and  Thermodynamics  of  Compressible  Fluid  Flow, 
Rowland  Press,  New  York,  1958. 
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In  considering  the  boundary  condition  problem,  guidance  car.  re  c-btained 
from  a  characteristics  analysis  of  the  inviscid  equations.  The  cr.arac  ttrri  s  t  i  cs 
analysis  reveals  how  boundary  data  is  propagated  into  the  solution  domain  by 
the  differential  equations  and  thus  can  be  used  to  determine  how  boundary  condi¬ 
tions  should  be  chosen  for  an  inviscid  flow  solution  procedure.  The  number  of 
boundary  conditions  required  at  any  boundary  point  is  equal  to  the  number  of 
characteristics  originating  at  the  point  which  transfer  information  from  the 
boundary  point  into  the  interior  flow  field.  Thus,  the  charac terist ics  analysis 
determines  the  number  of  boundary  conditions  to  be  applied  to  each  boundary 
point,  but  does  not  specify  particular  boundary  condition  equations.  The  invis¬ 
cid  characteristics  analysis  shows  that  two  boundary  conditions  are  required  at 
a  subsonic  inflow  boundary,  one  boundary  condition  is  required  at  a  subsonic 
outflow  boundary,  three  boundary  conditions  are  required  at  a  supersonic  inflow 
boundary  and  zero  boundary  conditions  are  required  at  a  supersonic  outflow 
boundary . 

Although  these  boundary  conditions  would  suffice  for  a  characteristics 
solution,  the  situation  becomes  more  complicated  when  a  finite  difference  method 
is  used  to  solve  the  equation  set.  If,  as  in  the  present  calculation,  central 
spatial  differences  are  used  to  represent  spatial  derivatives  and  the  equations 
are  written  at  each  interior  point,  then  three  conditions  must  be  written  at 
each  boundary  point  to  complete  the  resulting  set  of  algebraic  equations.  As 
indicated  by  the  characteristics  analysis,  the  differential  equations  allow 
specif ica t ion  of  three  conditions;  the  remaining  three  conditions  required  to 
close  the  set  of  equations  result  solely  from  the  use  of  central  differences 
and  they  are  termed  extraneous  conditions.  These  extraneous  boundary  conditions 
must  be  consistent  with  the  physical  flow  situation  or  spurious  numerical  solu¬ 
tions  may  result.  These  concepts  regarding  boundary  condition  specification 
are  tested  subsequently  by  solving  simple  problems.  In  the  test  calculations 
which  follow,  temporal  derivatives  are  represented  by  first-order  backward  time 
differences  and  spatial  derivatives  are  represented  by  second-order  accurate 
central  differences.  Following  Briley  and  McDonald  (Ref.  14)  the  equations  are 
solved  by  a  Linearized  Block  Implicit  (LBI)  scheme  where  nonlinear  terms  are 
linearized  by  Taylor  series  expansion  about  the  solution  at  a  known  time  level 
and  the  equations  are  solved  as  a  coupled  set.  Details  of  the  procedure  are 
given  in  Ref.  14. 


The  boundary  condition  test  case  for  use  v:ith  central  d:ff.* rtnevs 

is  first  explored  for  one-dinensional  subsonic  and  supersonic  uniform.  flows  at 
a  Reynolds  number,  Re^,  of  approximately  2.2  x  10^.  The  quantity  Re.  is  the 
Revnolds  number  of  the  problem  defined  as  Re „  =  oui/y  where  p,  u,  u  are  free 

X 

stream  density,  velocity  and  viscosity  respectively  and  Z  is  the  domain  length. 
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The  cell  Revnolds  number,  Re,  ,  in  these  cases  was  of  the  order  10  where 
Re  =  puax/p  and  Lx  is  the  grid  spacing.  Central  difference  solutions  may 
develop  spatial  oscillations  at  high  cell  Reynolds  numbers  unless  artificial 
dissipation  is  added  to  the  equations  either  explicitly  (by  the  addition  of  a 
dissipation-like  term)  or  implicitly  (e.g.,  one-sided  differencing).  For 
certain  model  problems  representing  a  balance  of  convection  and  normal  diffusion 
(Ref.  22),  the  maximum  cell  Reynolds  number  which  can  be  obtained  without  the 
appearance  of  spatial  oscillations  is  2,  and  many  numerical  procedures  require 
some  artificial  dissipation  if  the  cell  Reynolds  number  for  the  problem  being 
solved  and  grid  being  used  is  above  this  value.  This  point  will  be  considered 
in  some  detail  subsequently  when  artificial  dissipation  is  discussed.  The 
equations  solved  were  the  continuity  equation,  the  momentum  equation  and  the 
energy  equation  with  q  =  0.  A  nonuniform  distribution  of  the  dependent  vari¬ 
ables  was  specified  at  time  t  =  0  and  the  solution  marched  in  time;  the  exact 
solution  at  large  times  is  uniform  flow  with  constant  p,  u,  T°. 

As  discussed  previously,  the  characteristics  analysis  shows  that  for  sub¬ 
sonic  flow  two  conditions  must  be  specified  on  the  inflow  boundary  and  one  on 
the  outflow  boundary.  With  this  as  a  guide,  total  temperature  and  total  pres¬ 
sure  were  specified  on  the  inflow  boundary  and  static  pressure  was  specified 
on  the  outflow  boundary.  Since  spatial  derivatives  are  represented  by  central 
differences,  three  extraneous  boundary  conditions  are  required.  For  these  cases 
the  second  derivative  of  density  was  set  to  zero  on  the  inflow  boundary,  and 
second  derivatives  of  velocity  and  total  temperature  were  set  to  zero  on  the 

22.  Roache,  P.  J.:  Computational  Fluid  Dynamics.  Hermosa  Press, 

Albequerque,  1972. 
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dcwnstream  boundary.  This  set  of  boundary  conditions  is  referred 
subsonic.  It  should  be  noted  that  the  specified  extraneous  cone: 
consistent  with  the  exact  solution.  Additional  cases  were  run  wi 
conditions  chosen  for  supersonic  flow.  In  these  cases  total  temp 
velocity,  and  density  were  specified  at  the  inflow  boundary  and  t 
derivatives  of  these  quantities  were  set  to  zero  at  the  outflow  b 
set  of  conditions  is  termed  supersonic-supersonic.  Although  vise 
are  solved  in  this  test  case,  a  high  Reynolds  number  is  considers 
behavior  is  essentially  inviscid.  Therefore,  the  guidance  obtain 
inviscid  characteristics  analysis  would  be  expected  to  be  valid, 
cases  studied  are  given  in  TABLE  I 
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th  bour.  darv 
^rature , 
he  second 
oundary .  This 
ous  equations 
d  and  the  flow 
ed  from  the 
The  four 


TABLE  I  -  UNIFORM  FLOW  CASES 


Case 

M 

Boundary'  Conditions 

1 

<1 

Subsonic 

Subsonic 

2 

<1 

Supersonic  - 

Supersonic 

3 

>1 

Subsonic 

Subsonic 

4 

>1 

Supersonic  - 

Supersonic 

! 


The  test  cases  were  initiated  as  nonuniform  flows,  all  subsonic  or  all 
supersonic,  as  the  case  might  be,  which  satisfied  the  boundary  conditions,  and 
then  the  flow  was  allowed  to  develop  in  time.  Although  the  calculations  were 
run  at  a  high  cell  Reynolds  number  using  spatial  central  differences  with  no 
artificial  dissipation  term  added  for  numerical  stability,  the  entirely  subsonic 
flow  case  with  subsonic-subsonic  boundary  conditions  and  the  entirely  supersonic 
flow  case  with  supersonic-supersonic  boundary  conditions  converged  rapidly  to 
the  exact  uniform  solution.  However,  both  the  subsonic  flow  run  with  supersonic- 
supersonic  boundary  conditions  and  the  supersonic  flow  run  with  subsonic-subsonic 
boundary  conditions  developed  severe  spatial  oscillations  and  failed  to  converge. 
The  two  convergent  calculations  were  repeated  with  initial  conditions  specified 
which  did  not  satisfy  the  zero  derivative  boundary  conditions.  The  discrepancy 
between  the  initial  flow  and  the  specified  boundary  conditions  did  not  hinder 
convergence. 
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These  test  calculations  demonstrate  that  converged  solutions  can  be 
obtained  for  both  subsonic  and  supersonic  uniform  flow  using  central  -ratial 
differences  in  each  case,  provided  boundary  conditions  are  treated  adequately. 
One-sided  or  retarded  differences  were  not  required  to  obtain  a  converged 
supersonic  flow  solution.  These  boundary  conditions  used  in  the  converged 
calculations  are  consistent  with  an  inviscid  characteristics  analysis  and 
employ  exact  extraneous  boundary  conditions.  No  explicit  or  implicit  artificial 
dissipation  was  required  in  this  simple  case  of  unidirectional  flow  which  is 
uniform  in  the  steady  state. 

Nonuniform  Subsonic  Flow 

Further  test  calculations  were  made  for  one-dimensional  nonuniform  sub¬ 
sonic  flow.  As  in  the  subsonic  uniform  flow  case,  total  temperature  and  total 
pressure  were  set  as  upstream  boundary  conditions  and  the  density  gradient  at 
the  upstream  boundary  was  set  to  zero.  The  downstream  boundary  conditions 
consisted  of  specified  static  pressure  and  zero  derivatives  of  velocity  and 
total  temperature.  The  flow  was  accelerated  by  sped f ica t ions  of  a  heat  source 
distribution  in  the  energy  equation,  Eq .  (10);  however,  the  heat  sources  at 
both  ends  of  the  domain  were  brought  smoothly  to  zero  so  as  to  be  consistent 
with  the  zero  gradient  boundary  conditions  used.  A  dependent  variable  distri¬ 
bution  was  assumed  for  an  initial  flow  field  and  the  solution  allowed  to  develop 
in  time.  With  the  specified  heat  distribution,  the  solution  converged  and  upon 
convergence  the  values  of  the  dependent  variables  at  the  upstream  boundary  were 

W,  =  1.1062  W  r 
1  ref 

P,  «  .9898  P  . 

1  ref 

T°  =  1.0166  T°  , 

1  ref 

where  W  is  axial  velocity,  P  is  pressure  and  T°  is  total  temperature.  The  exit 
conditions  for  an  inviscid  Rayleigh  solution  with  the  above  inlet  conditions 
and  the  same  specified  heat  source  distribution  are  compared  to  the  computed 
Re  -  10^  solution  in  Table  II. 
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TABLE  II  -  NONUNIFORM  SUE  SONIC  FLOkT 

Rayleigh  solution  Present  solution 


W2 

1 . 2004 

1.2000 

P2 

.9771 

.9774 

T° 

2 

1.0919 

1.0938 

As  can  be  seen,  the  numerical  solution  at  high  Reynolds  number  is  in  excellent 
agreement  with  the  inviscid  Rayleigh  solution,  and  despite  the  very  high  cell 
Reynolds  number  convergence  was  obtained  readily  without  significant  spatial 

oscillations  using  central  spatial  differences  and  with  no  artificial  dissipa¬ 
tion  . 


Mixed  Subsonic-Supersonic  Flow 

The  next  flow  situation  investigated  was  characterized  by  positive  heat 
sources  in  the  upstream  region  of  the  flow  and  negative  heat  sources  (heat 
sinks)  in  the  downstream  region.  A  continuous  initial  flow  condition  which 
contained  subsonic  flow  in  the  region  of  the  heat  sources  and  supersonic  flow 
in  the  region  of  the  heat  sinks  was  chosen.  Upstream  boundary  conditions  of 
specified  total  pressure,  specified  total  temperature  and  zero  density  deriva¬ 
tive  were  set;  as  determined  earlier  these  conditions  appear  adequate  for  a 
subsonic  inflow  boundary.  The  downstream  boundary  conditions  specified  zero 
gradient  for  all  three  dependent  variables,  velocity,  density  and  total  tempera¬ 
ture;  as  determined  earlier  these  conditions  are  consistent  with  the  exact 
solution*  Thus,  for  the  mixed  flow  problem  both  inflow  and  outflow  boundary 
conditions  were  set  consistent  with  an  inviscid  characteristics  analysis  with 
extraneous  conditions  set  to  be  consistent  with  the  exact  solution.  The 
Reynolds  number  was  set  at  10^,  and  the  solution  was  then  allowed  to  develop 
with  time. 
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The  solution  first  approached  a  flow  field  consistent  with  an  inviscid 
Rayleigh  solution  having  a  throat  at  the  location  where  the  'neat  sources 
changed  from  positive  to  negative  values.  During  this  portion  of  the  calcula¬ 
tion,  the  change  in  the  dependent  variables  across  a  time  step  became  less 
than  one  percent,  with  the  time  step  being  the  time  required  for  a  particle 
traveling  at  sonic  velocity  to  pass  from  one  grid  point  to  another.  However, 
the  calculation  was  allowed  to  continue  and  the  solution  began  to  diverge  with 
an  "expansion  shock"  appearing  across  the  sonic  region.  The  difference  equa¬ 
tions  across  the  shock  satisfied  the  momentum  and  continuity  equations  with  a 
T°  change  consistent  with  the  heat  source.  Thus,  the  solution  shifted  abruptly 
from  the  subsonic  branch  of  the  Rayleigh  curve  to  the  supersonic  branch  of  the 
Rayleigh  curve.  This  solution  represents  a  physically  unrealizable  flow  field 
(e.g.  Ref.  21).  As  the  calculation  continued,  the  strength  of  the  expansion 
shock  increased  with  the  subsonic  region  being  driven  towards  stagnation  condi¬ 
tions  and  the  supersonic  region  being  driven  toward  very  high  Mach  numbers. 
Convergence  of  this  solution  to  a  steady  state  did  not  appear  possible  and  the 
solution  was  discontinued.  It  appeared  that  without  something  akin  to  artifi¬ 
cial  dissipation  a  solution  could  not  be  obtained  at  this  Reynolds  number  with 
the  boundary  conditions  taken  from  the  inviscid  characteristics  analysis. 

Test  Calculations  with  Artificial  Dissipation 
Second  Order  Dissipation 

The  "expansion  shock"  discussed  in  the  previous  section  is  likely  to 
occur  only  at  high  Reynolds  numbers.  If  large  viscous  terms  were  present, 
they  would  tend  to  dissipate  the  discontinuous  flow  region.  Two  general 
methods  exist  for  adding  nonphysical  dissipative  mechanisms  to  the  equations. 

In  the  first  method,  the  truncation  error  associated  with  the  spatial  difference 
representation  represents  a  dissipative  type  term.  An  example  of  this  approach 
is  the  retarded  difference  representation  (e.g.  Ref.  16).  In  the  second  general 
method,  dissipation  is  added  explicitly  through  an  additional  term.  Since  the 
present  approach  utilizes  spatial  central  differences  which  do  not  have  a 
dissipative  type  truncation  error  associated  with  them,  the  present  approach 
utilizes  an  explicit  artificial  dissipation  whose  magnitude  can  be  precisely 
controlled. 
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Several  different  methods  cf  adding  artificial  dissipative  terns  Lave 
been  proposed,  and  all  atter.pt  to  add  dissipative  phenomena  only  in  regions 
of  the  flow  where  the  artificial  terns  are  required  either  for  numerical  sta¬ 
bility  or  to  suppress  the  short  wavelength  oscillations  generated  by  the  dis¬ 
continuity.  Three  commonly  used  artificial  dissipation  devices  are  second- 
order  dissipation  (artificial  viscosity),  fourth-order  dissipation  and 

dissipation  coefficient  based  upon  pressure  derivative.  In  second-order 

2  ? 

dissipation,  a  second-order  term  of  the  form  q?  is  added  to  the 

governing  equation  where  a  is  a  dissipation  coefficient  and  <*  is  a  flow  var¬ 
iable.  As  usually  applied,  fourth-order  dissipation  adds  an  entirely  new  type 

4  4 

of  term  to  the  governing  equations,  of  the  form  -  a  o  4>/3  x  .  In  the  final 
type  of  artificial  dissipation  considered,  the  additional  dissipative  term  is 
made  proportional  to  the  second  derivative  of  pressure.  In  all  cases,  the 
question  arises  as  to  whether  the  addition  of  the  artificial  dissipation  terms 
would  destroy  the  essential  features  of  solutions  which  contain  very  rapid  un¬ 
resolved  spatial  changes  in  the  dependent  variables.  Therefore,  a  study  of  the 
effect  of  adding  artificial  dissipation  to  the  governing  set  of  equations  was 
made. 

In  pursuing  this  study  and  considering  the  results  which  will  be  presented, 
several  items  should  be  kept  in  mind.  First  of  all,  the  one-dimensional  problem 
which  serves  as  the  test  problem  represents  a  balance  between  convective  and 
streamwise  diffusive  phenomena.  These  are  the  important  phenomena  in  shock  flow 
considerations,  but  these  are  not  the  dominant  phenomena  in  other  important 
applications  such  as  shock-free  regions  of  shear  flow  where  transverse  diffusion 
may  be  critical.  Secondly,  the  equations  are  being  solved  by  an  LBI  method 
(Ref.  14)  with  second-order  accurate  central  differencing  used  throughout. 

Thus,  (except  for  the  negligible  axial  diffusion  due  to  viscous  terms  in  this 
10^  Reynolds  number  flow)  the  artificial  viscous  terms  represent  the  only  dis¬ 
sipative  terms  present.  Other  methods  may  contain  more  or  less  spatial 
dissipation  due  to  truncation  error  which  could  be  of  higher  or  lower  order  and, 
therefore,  other  solution  procedures  might  behave  differently.  Finally,  the 
present  analysis  utilizes  density,  velocity  and  total  temperature  as  dependent 
variables.  Other  procedures  may  choose  alternate  sets  of  dependent  variables, 
with  consequent  changes  in  the  observed  effect  of  the  artificial  dissipation. 
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In  the  second-order  dissipation  formulation  a  value  is  chosen  for  the  dissipation 
parameter,  o,  and  this  choice  leads  to  a  field  specification  of  and  with 
the  stipulation  that  neither  be  negative.  High  values  of  o  imply  relatively 
large  artificial  dissipation  terms  and  low  values  imply  relatively  small 
artificial  dissipation  terms. 

The  subsonic- inf low  supersonic-outflow  case  described  previously  then  was 
recalculated  with  the  artificial  dissipation  terms  included  for  a  series  of 
dissipation  parameters.  As  shown  in  Fig.  1,  the  calculation  made  with  a  dis¬ 
sipation  parameter  value  of  .5  was  severely  damped  and  differed  significantly 
from  the  inviscid  Rayleigh  solution.  It  is  believed  that  in  view  of  the  physical 
Reynolds  number  of  10^  there  should  be  little  difference  between  the  Rayleigh 
solution  and  an  accurate  numerical  solution.  When  an  o  -  .025  condition  was 
imposed,  the  calculated  solution  agreed  very  well  with  the  inviscid  solution. 

The  last  case  considered  set  o  at  a  value  of  .0025  and  although  convergence 
was  obtained,  a  physically  unrealistic  "expansion  shock”  similar  to  that 
obtained  earlier  without  any  artificial  dissipation  terms  present  was  predicted 
at  the  throat.  Thus,  the  results  of  Fig.  1  inidcate  that,  for  the  simple  case 
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presented,  an  artificial  viscosity  formulation  can  give  physically  r  a i  c 
results  with  a  suitable  :*oice  of  dissipation  parameter.  Vc.en  the  c  :  -  l : :  a:  i .  n 
parameter  is  set  to  0.5,  the  results  are  smeared  so  that  they  significantly 
differ  from  the  desired  inviscid  solution.  However,  at  lower  values  of  the 
dissipation  parameter,  0  -  .025,  a  solution  which  shows  good  agreement  with  the 
known  inviscid  solution  is  obtained  and  finally,  at  the  lowest  dissipation 
parameter  investigated,  a  =  .0025,  a  physically  unrealistic  discon t inui tv 
appears  in  the  solution.  When  considering  the  results  of  this  test  calculation, 
it  should  be  emphasized  that  the  model  problem  is  dominated  by  a  balance  between 
convective,  pr assure  and  streanw'ise  diffusion  terms.  It  does  not  address  the 
shock-free  shear  layer  problem  where  the  balance  is  dominated  by  convective, 
pressure  and  transverse  diffusion  terms. 

A  different  flow  situation  can  be  obtained  with  the  same  distribution  of 
heat  sources  and  sinks  by  imposing  a  downstream  static  pressure  boundary 
condition.  However,  once  the  total  temperature  and  static  pressure  at  the 
sonic  throat  are  determined,  and  if  the  heat  source  distribution  downstream  of 
the  throat  is  specified,  then  the  Rayleigh  solution  admits  only  two  possible 
exit  pressures.  The  first  is  the  supersonic  pressure  shown  in  Fig.  1;  the 
second  is  a  subsonic  pressure  which  occurs  when  a  shock  is  present  at  any 
location  between  the  throat  and  the  downstream  boundary.  Unlike  the  case  of 
flow  in  a  nozzle  where  the  nozzle  exit  pressure  varies  with  the  shock  position, 
the  Rayleigh  line  solution  gives  the  same  downstream  pressure  for  all  shock 
locations . 

For  the  conditions  of  the  model  problem  the  permissible  downstream  subsonic 
pressure  is  1.19  times  the  static  pressure  at  the  throat,  p*  (Ref.  21).  A  dis¬ 
sipation  parameter  study  was  made  for  one-dimensional  flow  with  the  pressure  at 
the  downstream  boundary  being  1.19p*,  and  the  results  are  shown  in  Fig.  2. 

The  calculation  made  for  0  -  .025  predicted  a  shock  to  occur  when  the  flow  en¬ 
countered  an  adverse  pressure  gradient.  The  shock  was  spread  over  three  grid 
points  and  except  for  a  small  post-shock  overshoot ,  the  result"  compared  favor¬ 
ably  with  a  Rayleigh  line  calculation  that  contains  a  shock  at  the  same  location. 
The  c  =  .025  calculation  converged  over  most  of  the  flow  field;  however,  at  the 
end  of  the  run  temporal  changes  (oscillations)  were  still  occurring  in  the  im¬ 
mediate  vicinity  of  the  shock.  The  calculation  with  0  =  .5  once  again  converged 
but  showed  a  considerable  discrepancy  with  the  inviscid  solution  as  a  result  of 
excess  dissipation.  In  fact,  the  0  -  .5  solution  did  not  even  contain  a  super- 
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sz:\ic  Finally,  liie  j  =  .0025  calculation  became  unstable  afc  3  la  roc 

:r*.  shock  ovo  reboot  occurred  and  the  cal  dilation  w.>uld  not  ccr.vcrpc. 

Since  only  two  downstream  static  pressures  are  permissible  when  the  irvis- 
cid  flow  is  sonic  at  the  throat,  it  is  instructive  to  examine  the  results  cf 
calculations  made  at  high  Reynolds  numbers  when  downstream  boundary  conditions 
inconsistent  with  the  inviscid  flow  are  specified.  Figure  3  shows  the  results 
of  calculations  made  when  the  downstream  pressure  is  set  at  O.S12p*;  this  down¬ 
stream  pressure  is  lower  than  the  permissible  inviscid  supersonic  value.  The 
c  =  .5  calculation  proceeds  smoothly  into  this  boundary  condition  through  a 
sudden  acceleration  in  the  downstream  portion  of  the  calculation  region.  The 
3  =  .023  calculation  shows  severe  spatial  oscillations  and  although  the  calcu¬ 
lation  does  not  appear  to  diverge,  the  calculation  would  not  converge  to  a 
steady  state  in  the  allotted  time.  The  results  of  a  calculation  made  for 
-  P*  are  shown  in  Fig.  4.  This  value  of  back  pressure  is  between  permis¬ 
sible  subsonic  and  supersonic  values  and  again  represents  a  physically  unreal¬ 
izable  inviscid  flow.  Once  again  the  a  =  .5  calculation  appears  smooth  and 
satisfies  the  boundary  condition.  However,  the  o  =  .025  calculation  shows 
severe  spatial  oscillations  and  although  once  again  it  did  not  diverge,  it  did 
not  achieve  steady  state.  Thus,  for  this  simple  flow  the  o  =  .5  criterion 
accepts  the  imposition  of  physically  unrealistic  inviscid  boundary  conditions 
without  exhibiting  any  obvious  inconsistent  behavior. 

Fourth-Order  Dissipation 

The  second-order  dissipation  approach  represents  one  technique  for  adding 

artificial  dissipation  to  the  system  of  viscous  flow  equations;  however,  it 

represents  only  one  of  several  methods  which  have  been  suggested.  One  alterna- 

4  4 

tive  method  is  the  fourth-order  dissipation  approach  where  the  term  (-0^3  <p/2x  ) 
is  added  to  the  right-hand  sides  of  the  equations;  <f>  represents  p,  u,  T°  for 
the  continuity,  momentum  and  energy  equations,  respectively.  Such  an  approach 
has  been  used  by  several  investigators  (e.g.  Ref.  23)  sometimes  in  conjunction 


23.  Beam,  R.  M.  and  R.  F.  Warming:  An  Implicit  Factored  Scheme  for  the 

Compressible  Navier-Stokes  Eqs.  Proceedings  of  AIAA  Third  Computational 
Fluid  Mechanics  Conference,  1977. 
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v  i  :  h  second-order  dissipation.  In  sor.t  nrpr^eches  t  coeffic’ont  -  i?  a 
function  of  the  time  step  (e.g.,  Ref.  2  3);  however,  ::  .  is  a  fwnciion  if  t:.v 

time  step,  the  steady  state  finally  reached  will  be  dependent  upon  the  tine 
step  used  in  the  calculation;  therefore,  was  made  independent  of  time  step 
in  the  present  inves t igat ion .  There  are  two  general  methods  of  including 
fourth-order  dissipation  in  the  numerical  procedure.  In  the  implicit  treat¬ 
ment,  the  size  of  the  matrix  being  inverted  increases  from  a  block  tridiagonal 
to  a  block  pent ^diagonal  form  and  the  cost  of  matrix  inversion  mere  than 
doubles.  This  large  increase  in  the  matrix  inversion  time  is  difficult  to 
justify  on  the  basis  of  artificial  dissipation  alone  and,  therefore,  the 
implicit  inclusion  of  fourth-order  dissipation  was  not  considered  further. 

The  second  method  considered  includes  the  additional  terms  in  the  difference 
equations  in  an  explicit  manner.  Since  this  does  not  increase  the  computer 
run  time  by  any  significant  amount,  this  explicit  approach  was  evaluated. 

The  one-dimensional  shock  calculation  was  made  for  three  different  values  of 
a,  ,  a,  =  0.03,  0.20,  1.0.  In  all  cases  the  calculation  included  an  insignifi- 
cant  amount  of  second-order  dissipation  with  o  held  at  .00025.  When  the 
smallest  value  of  was  used,  cl^  =  0.03,  the  resulting  calculation  did  not 
converge,  and  contained  both  a  large  "expansion  shock"  near  the  throat 
and  a  large  compression  shock,  near  the  downstream  boundary.  The  calculation 
made  using  the  large  value  of  *  1,  also  gave  very  poor  results  as  it  devel¬ 
oped  large  diverging  spatial  and  temporal  oscillations.  This  adverse  behavior 
shown  by  the  large  calculation  may  be  related  to  violation  of  an  effective 
explicit  stability  limit  associated  with  the  explicit  representation  of  the 
dissipation  terra. 

The  final  fourth-order  dissipation  case  considered  was  run  with  =  0.20. 
The  results  of  this  calculation  are  shown  in  Fig.  5.  As  can  be  seen  in  Fig.  5, 
the  results  are  not  very  different  from  the  second-order  o  *  .025  case,  the 
solution  was  unchanging  over  most  of  the  flow  field  but  at  the  end  of  the  run 
some  temporal  changes  were  still  occurring  in  the  vicinity  of  the  shock. 

Pressure  Derivative  Damping  Coefficient 

The  final  approach  considered  follows  the  'shock  smearing1  technique  used 

3  ’  2  2  2  2 

in  Ref.  24  and  adds  the  term  a, (Ax)  (|u|+c)/4p  3  p/Bx  B  $/3x  to  the  right-hand 

^  3  2  2 

side  of  each  equation.  The  factor  a^(Ax)  (|u|+c)/4p  B  p/Rx  may  be  interpreted 
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2  S  3 


1 


/iscous  diffusion  coefficient  whose  strength  depends  upor.  the  h'w: 


sure  crauie: 


Viewed  in  this  manner  tht 


ior.al  tern,  r.av  b< 


a  cuusil inear  second-order  tern,  which  can  be  treated  implc fitly.  Two  calcu¬ 
lations  were  run  with  this  formulation;  the  first  assumed  a,  =  0.0005  and  the 

4 

second  set  a,  at  0.02.  In  the  a.  -  0.005  calculation,  steady  state  was  never 

obtained,  as  the  calculation  developed  strong  spatial  oscillations  which  varied 

with  time.  The  second  calculation  made  for  a.  =  0.02  appeared  much  smoother; 

4 

however,  the  predicted  shock  was  smeared  considerably  more  than  the  a  -  .025 
=>econd-orcer  dissiDation  calculation. 


Mass  Flux  Damping 


The  test  cases  presented  in  Figs.  1-6  all  consider  the  dependent  variable 
to  be  p,  u,  T° ,  however,  another  possibility  is  to  choose  p,  pu,  and  T°  as 
dependent  variables.  In  this  latter  case  the  artificial  dissipation  terms  are 
constructed  so  as  to  suppress  the  short  wavelength  spatial  oscillations  in  the 
quantities  p,  pu,  and  T°.  For  the  case  of  second  order  dissipation  this 
consists  of  adding  terms 


3  p  a  pu 


d2  32T° 


3x 


to  the  governing  equations  with  d^  and  d2  being  chosen  so  as  to  satisfy  a  dis¬ 
sipation  parameter  criterion.  This  set  of  artificial  dissipation  quantities  is 
termed  SET  B  whereas  the  original  set  based  upon  the  terms 


1  9  p 

P  a  2 
9x 


d  3ju 
3x2 


,  2  o 

d  9  T 

9x 


is  termed  SET  A.  Results  of  the  ’SET  B’  artificial  dissipation  formulation 
along  with  a  Rayleigh  line  solution  are  presented  in  Fig.  7  for  three  values  of 
dissipation  parameter.  As  shown  in  Fig.  7,  the  solutions  for  o  of  .5  and  .1 
both  contain  excessive  damping  as  neither  properly  represents  the  shock  region. 
The  results  for  a  «  .05  shows  better  agreement  with  the  Rayleigh  solution, 
however,  the  computed  solution  contains  relatively  strong  spatial  oscillations. 
A  comparison  of  ’SET  A’  and  ’SET  B’  damping  formulations  for  dissipation 
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parameters  of  .05  and  .025  are  shewn  in  Figs.  S  and  9.  As  c^n  be  re^n  :r jr  t 
figures,  a  better  numerical  representation  of  the  problem  apt-ear?  u-  result 
from  the  ’SET  A*  artificial  dissipation  formulation  (p,  u,  1°)  at  c  =  .025. 

Summary 

This  one-dimensional  nodel  problem  has  served  to  focus  upon  properties 
of  artificial  dissipation  techniques  which  are  relevant  to  the  transonic 
shock  flow  field.  The  problem  itself  represents  a  case  where  the  flow  field 
equations  are  balanced  by  convective-pressure  and  axial  diffusion  phenomena. 

As  in  the  transonic  cascade,  the  shock  results  from  the  specif ica t ion  of  a 
downstream  static  pressure  condition  and,  therefore,  the  mechanism  for  shock 
generation  in  the  transonic  cascade  problem  of  interest  and  in  the  one¬ 
dimensional  model  problem  are  very  similar.  With  these  considerations  in  mind, 
results  of  the  model  problem  analysis  can  give  guidance  to  the  transonic 
cascade  shock  capturing  problem. 

The  numerical  results  for  the  model  problem  can  be  assessed  by  comparison 
to  one-dimensional  inviscid  theory  for  flow  with  heat  sources.  Although  this 
Rayleigh  solution  is  not  unique  in  terms  of  the  shock  position,  once  this 
position  is  specified,  the  analytic  solution  can  be  used  to  evaluate  the  artifi¬ 
cial  dissipation  models.  For  (i)  the  solution  procedure  used  and  (ii)  the 
problem  considered,  second-order  dissipation  with  a  dissipation  parameter,  a, 
equal  to  0.025,  gives  very  good  agreement  with  the  ’shocked  Rayleigh  solution1. 
Obviously,  extrapolat ion  of  these  results  to  the  complicated  two-dimensional 
transonic  cascade  problem  must  be  somewhat  speculative.  However,  since  the 
shock  in  the  transonic  case  is  aligned  fairly  closely  with  a  coordinate  line, 
it  is  reasonable  to  expect  that  an  accurate  shock  capturing  procedure  could  be 
obtained  by  adding  second-order  artificial  dissipation  to  the  governing  equa¬ 
tions  with  the  dissipation  parameter  taken  to  be  the  order  of  0.025. 


27 


III. 


CASCADE  ANALYSIS 
Coordinate  System 

With  increasing  development  of  computational  techniques  many  of  the  pro¬ 
cedures  originally  developed  to  solve  viscous  flow  equations  in  simple  coordi¬ 
nate  systems  are  being  applied  to  realistic  geometric  configurations.  Such  ap¬ 
plications  require  the  development  of  suitable  computational  coordinate  systems. 
For  a  coordinate  system  to  be  computationally  suitable  it  must  satisfy  condi¬ 
tions  of  (i)  generality,  (ii)  smoothness,  (iii)  resolvability  and  (iv)  allow 
easy  implementation  of  boundary  conditions.  Considering  these  conditions  one 
by  one,  a  coordinate  generation  technique  should  be  fairly  general  allowing 
specification  of  airfoil  shapes  and  cascade  configurations  of  practical  in¬ 
terest.  Secondly,  the  technique  must  generate  coordinate  systems  sufficiently 
smooth  to  allow  reliable,  accurate  calculations.  This  requirement  is  more 
stringent  than  simply  requiring  the  existence  of  a  given  number  of  coordinate 
transformation  derivatives;  the  coordinates  must  not  have  large  changes  in  the 
coordinate  geometry  data  between  grid  points  for,  if  large  changes  exist,  the 
solution  accuracy  will  deteriorate.  An  assessment  of  the  smoothness  requires 
a  detailed  examination  of  the  metric  data.  Concerning  the  third  requirement, 
the  grid  generation  technique  must  allow  high  grid  resolution  in  specified 
regions  of  the  flow  field.  For  example,  the  wall  region  of  viscous  flows  is 
characterized  by  rapid  changes  in  flow  variables.  Therefore,  the  computational 
grid  should  have  a  large  number  of  points  in  the  viscous  wall  region.  Other 
regions  in  the  flow  are  regions  of  slow  changes  in  flow  variables  and  these 
regions  should  have  a  relatively  sparse  computational  grid  associated  with  them. 
Finally,  the  technique  must  allow  accurate  specif icat ion  of  the  required  boundary’ 
conditions.  If  flow  on  a  given  line  is  to  satisfy  a  no-slip  condition,  then  it 
is  important  that  this  no-slip  line  be  a  coordinate  line.  The  specification  of 
boundary  conditions  on  a  boundary  line  which  is  not  a  coordinate  line  may  pre¬ 
sent  a  difficult  problem  for  Navier-S tokes  analyses.  Although  finite  difference 
molecules  can  be  constructed  for  boundary  points  which  do  not  fall  upon  coordi¬ 
nate  lines,  such  molecules  may  have  an  unacceptable  truncation  error  associated 
with  them.  Further,  it  should  be  noted  that  the  boundary  condition  problem 
arises  in  both  viscous  and  inviscid  analyses;  however,  the  problem  is  consid¬ 
erably  more  severe  in  viscous  flows  where  no-slip  conditions  on  solid  walls 
can  combine  with  boundary  condition  truncation  error  to  produce  numerical 
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solutions  which  are  both  qualitatively  and  quantitatively  in  error.  Similarly, 
if  a  given  line  is  a  periodic  line,  it  is  also  important  that  tine  periodic 
line  be  a  coordinate  line  since  boundary  conditions  are  applied  here.  In 
addition,  application  of  the  periodic  condition  is  greatly  facilitated  if  the 
coordinate  metric  data  is  continuous  across  the  periodic  line.  If  the  generated 
metric  coefficients  are  d iscontinuous  across  periodic  lines,  a  device  such  as 
one-sided  differencing  would  be  required  to  represent  the  equations  at  this 
location  with  a  consequent  loss  of  accuracy. 

Among  the  coordinate  system  generation  techniques  currently  used  are  con¬ 
formal  mapping  techniques,  techniques  based  upon  a  solution  of  Poisson’s  equa¬ 
tion  and  constructive  techniques.  An  example  of  a  conformal  mapping  technique 
applied  to  the  cascade  problem  is  given  by  Ives  and  Liutermoza  (Ref.  3). 

Examples  of  coordinate  systems  generated  by  the  solution  of  a  Poisson  equation 
can  be  found  in  the  work  of  Thompson  and  his  coworkers  (e.g.  Ref.  25)  or 
Steger  and  Sorenson  (Ref.  26).  The  third  approach  to  the  coordinate  system 
problem  is  the  constructive  approach  used  by  Gibeling,  Shamroth  and  Eiseman 
(Ref.  27)  for  the  isolated  airfoil  problem  and  by  Eiseman  (Refs.  28  and  29) 
for  the  cascade  problem.  Application  to  the  isolated  airfoil  problem  is 
discussed  by  Shamroth  and  Levy  (Ref.  30)  and  a  description  of  its  application 
to  the  cascade  problem  is  discussed  by  Shamroth,  Gibeling  and  McDonald  (Ref.  13). 


25.  Thompson,  J.  F.,  F.  C.  Thames  and  C.  W.  Mastin:  Boundary  Fitted  Curvi¬ 
linear  Coordinate  Systems  for  Solution  of  Partial  Differential  Equations 
on  Fields  Containing  Any  Number  of  Arbitrary  Two-Dinensional  Bodies. 

NASA  CR-2729 ,  July  1977. 

26.  Steger,  J.  L.  and  R.  L.  Sorenson:  Automatic  Mesh  Point  Clustering  Near 

a  Boundary  in  Grid  Generation  with  Elliptic  Partial  Differential  Equations. 
Journal  of  Computational  Physics,  Vol.  33,  1979,  p.  405. 

27.  Gibeling,  H.  J.  ,  S.  J.  Shamroth  and  P.  R.  Eiseman:  Analysis  of  Strong- 
Interaction  Dynamic  Stall  for  Laminar  Flow  on  Airfoils.  NASA  CR-2969, 

April  1978. 

28.  Eiseman,  P.  R. :  A  Coordinate  System  for  a  Viscous  Transonic  Cascade 
Analysis.  United  Technologies  Research  Center  Report  R76-912149-4 ,  1976. 

29.  Eiseman,  P.  R. :  A  Coordinate  System  for  a  Viscous  Transonic  Cascade 
Analysis.  J.  Comp.  Physics,  Vol.  26,  March  1978,  pp.  307-338. 

30.  Shamroth,  S.  J.  and  R.  Levy:  Description  of  a  Constructive  Coordinate 
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Although  the  analysis  described  in  Ref.  13  in  principle  could  be  used  for 


cascade  of  arbitrary  airfoils,  the  coordinate  generation  code  vas  liritL'i  to 
unst3ggered  cascades  of  symmetric,  uncambered  airfoils.  The  cede  has  been 


generalized  by  Kim  and  Shamroth  (Ref.  31)  to  allow  specification  of  more  arbi¬ 


trary  airfoil  cascades  and  it  is  this  generalized  analvsis  which  was  used  in 


the  present  effort. 

In  brief,  the  constructed  coordinates  are  still  obtained  from  the  four 
loop  system  described  in  Ref.  13.  However,  with  the  generalization  described 
in  Ref.  31,  the  outer  loop  follows  the  camber  line  over  the  covered  portion  of 
the  cascade  and  the  blade  surface  itself  is  input  via  discrete  data  points. 

In  addition,  the  coordinate  orthogonality  condition  at  the  airfoil  surface  has 
been  relieved  in  the  vicinity  of  the  airfoil  leading  edge  and  in  this  region 
the  strict  orthogonality  condition  has  been  replaced  by  the  requirement  of  a 
smooth  variation  in  the  coordinate  angle. 

A  plot  of  the  coordinate  system  used  in  the  cascade  calculations  of  the 
present  report  is  shown  in  Fig.  10  where  the  calculation  utilized  a  30  x  113 
grid.  The  actual  inner  curve  representing  the  airfoil  is  a  smooth  curve 
defined  by  300  specif icat ion  points.  Only  selected  coordinate  lines  are  shown 
on  this  plot  with  pseudo-radial  coordinate  lines  in  the  vicinity  of  the  leading 
edge  and  pseudo-azimuthal  lines  in  the  vicinity  of  the  airfoil  leading  edge 
being  omitted  for  clarity.  To  resolve  the  turbulent  boundary  layer  viscous 
sublayer  grid  pseudo-radial  spacing  in  the  immediate  vicinity  of  the  blade  was 
.00004  chords.  The  pseudo-az imuthal  spacing  in  the  vicinity  of  the  blade  lead¬ 
ing  edge  is  0.0046  chords,  as  resolution  of  the  suction  peak  was  believed  to 
be  a  critical  factor  for  the  streamvise  mesh.  In  addition,  the  grid  extends 
until  two  chords  downstream  of  the  blade  trailing  edge. 


The  Governing  Equations 


A  sufficient  set  of  governing  equations  to  solve  the  viscous  cascade 
problem  for  the  flow  conditions  of  interest  is  the  ensemble  averaged  Navier- 
Stoices  equations  along  w:  th  the  usual  continuity  equation  and  an  energy  equa¬ 
tion.  However,  there  still  remain  decisions  to  be  made  on  the  choice  of 


31.  Kim,  Y.  N.  and  S.  J.  Shamroth:  Revised  Coordinate  Generation  Program 
for  a  Cascade  of  Arbitrarily  Shaped  Airfoils.  SRA  Rpt.  81-1,  1981. 
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dependent  variables  and  the  precise  forn  c:  these  equations.  T::c  r-r.bler.  is 
further  complicated  by  the  presence  of  a  complex  geometry.  Mar.v  previous 
Navier-Scokes  analyses  have  concentrated  upon  flows  in  Cartesian  or  cvlindrical 
geometries  and  in  these  cases  the  form  of  equations  is  considerably  simpler 
than  the  form  which  occurs  in  the  general  case.  The  presence  of  complex 
geometries  has  required  a  reevaluation  of  the  choice  of  equation  form  and 
dependent  variables,  and  such  a  reevaluation  has  been  made  by  Shamroth  and 
Gibeling  (Ref.  13).  As  a  result  of  the  reevaluation.  Ref.  13  suggests  the 
governing  equations  be  transformed  from  the  Cartesian  coordinates  x,y  to  a 
new  set  of  coordinates  £ ,  r, ,  where 


£  *  £  (  x  ,  y ,  t ) 
V  ■  ,y,t) 


T  •  t 


and  then  written  as 
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This  formulation  is  termed  the  quasi-linear  forr.  and  has  bt-j-  v-l-c  - 

fully  for  a  number  of  cascade  and  airfoil  calculations  (Refs.  13,  31'  and  33) 
for  both  laminar  and  turbulent  flow  in  the  subsonic  regime,  however,  the  formu¬ 
lation  has  not  yet  been  used  in  the  transonic  regime.  For  transonic  flow  two 
additional  items  must  be  considered;  these  are  the  choice  of  dependent  variables 
and  the  use  of  the  "quasi-linear"  equation  form  Eq .  (14)  rather  than  the  strong 

ccnservat ion  form,  which  is  written  as 


dw /D  d 
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Considering  the  first  of  these  items,  when  transonic  flow  is  considered,  two 
obvious  choices  for  dependent  variable  sets  arise.  These  are  the  primitive 
set  of  dependent  variables,  p,  u,  v  and  the  mass  flux  set  of  dependent  var¬ 
iables  p,  pu,  pv.  Within  the  framework  of  the  LBI  procedure  (Ref.  14)  the 
steady  converged  solution  obtained  using  either  set  of  dependent  variables 
will  be  the  same;  however,  differences  will  occur  if  an  artificial  dissipa¬ 
tion  terra  is  added  to  the  equations.  Therefore,  for  convenience  the  same 
primitive  set  of  dependent  variables  is  used  here  as  was  used  in  previous 
subsonic  calculations  (Ref.  13).  Furthermore,  based  upon  the  results  of  the 
one-dimensional  Rayleigh  problem  with  heat  addition,  second-order  dissipation 
based  upon  the  variables  p,  u,  v  is  used.  The  second  item  to  be  considered 
concerns  the  choice  of  Eq .  (14)  or  Eq.  (16)  and  results  of  good  quality  have 

been  obtained  for  steady  flows  using  the  ’'quasi-linear”  formulation,  Eq .  (14) 

for  subsonic  airfoil  and  cascade  flows  (Refs.  13,  31  and  32).  Further,  since 
the  metric  coefficients  are  independent  of  the  flow  being  transonic  and  since 


32.  Shararoth,  S.  J.  and  H.  J.  Gibeling:  A  Compressible  Solution  of  the 
Navier-Stokes  Equations  for  Turbulent  Flow  abc  ut  an  Airfoil.  NASA 
CR-3183 ,  1979.  (Also  AIAA  Paper  79-1543). 

33.  Shamroth,  S.  J.  and  H.  J.  Gibeling:  Analysis  of  Turbulent  Flow  about 
an  Isolated  Airfoil  with  a  Time-Dependent  Navier-Stokes  Procedure. 
Paper  presented  at  AGARD  Specialists  Meeting  on  Boundary  Layer  Effects 
on  Unsteady  Airloads,  Aix-en-Provence,  France,  September  1980. 
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results  obtained  via  Eq .  (14)  are  less  sensitive  to  the  method  used  tc  ^va'uaie 
the  metric  coefficients  than  are  results  obtained  via  Eq.  (16),  Eq.  (K;,  the.- 
quasi-linear  form  was  used  in  the  present  effort. 

Boundary  Conditions 

The  final  item  to  be  presented  prior  to  performing  cascade  calculations 
concerns  boundary  conditions.  The  boundary  conditions  used  in  the  present 
calculations  followed  the  suggestion  of  Briley  and  McDonald  (Ref.  34);  these 
were  used  successfully  in  the  previous  cascade  study  described  in  Ref.  13. 

For  the  cascade  system  shown  in  Fig.  10,  A3  and  CD  are  periodic  boundaries 
and  periodic  conditions  are  set  here.  It  should  be  noted  that  specification 
of  periodic  boundary  conditions  is  simplified  considerably  if  the  metric  data 
is  continuous  at  periodic  lines.  If  the  metric  data  is  not  continuous,  devices 
such  as  special  difference  molecules  or  one-sided  differencing  is  required  and 
in  general  these  devices  introduce  logical  complexity  and  may  cause  the  solu¬ 
tion  to  suffer  a  loss  of  accuracy.  The  present  grid  is  orthogonal  at  the 
periodic  line  and,  therefore,  the  metric  coefficients  are  continuous  across 
this  line.  With  the  present  grid  construction  process,  application  of  the 
periodic  boundary  conditions  is  a  straight-forward  task. 

Specifications  of  upstream  and  downstream  conditions  is  somewhat  more 
difficult.  For  an  isolated  cascade,  boundary  conditions  for  the  differential 
equations  may  be  known  at  both  upstream  infinity  and  downstream  infinity. 
However,  since  computation  economics  argues  for  placing  grid  points  in  the 
vicinity  of  the  cascade  and  minimizing  the  number  of  grid  points  far  from  the 
cascade,  the  upstream  and  downstream  computational  boundaries  should  be  set  as 
close  to  the  cascade  as  is  practical.  In  addition,  with  the  particular  body- 
fitted  coordinates  used  as  the  upstream  boundary  moves  further  upstream,  the 
angle  between  pseudo-radial  and  pseudo-azimuthal  coordinate  lines  becomes 
smaller.  Decreasing  the  coordinate  angle  causes  the  coordinate  system  to 
become  less  well-conditioned  and  increases  the  role  of  cross-derivative  terms 
in  the  equations.  Both  of  these  characteristics  could  be  detrimental  to  the 


34.  Briley,  V.1.  R.  and  H.  McDonald:  Computation  of  Three-Dimensional  Horse¬ 
shoe  Vortex  Flow  Using  the  Navier-Stokes  Equations.  7th  International 
Conference  on  Numerical  Methods  in  Fluid  Dynamics,  June  1980. 
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present  numerical  procedure  ar.d,  therefore,  they  also  :r:ue  :" :  r  placing  :he 
upstream  boundary  as  close  to  the  cascade  as  possible.  However,  v;h-:.  the 
upstream  boundary  is  placed  close  to  the  cascade,  most  flow  function  renditions 
on  the  boundary  will  not  be  known. 

In  the  present  approach  the  suggestion  of  Ref.  34  is  followed  which  sets 
total  pressure  on  boundary  BC.  Unless  boundary  BC  is  very  far  upstream,  the 
flow  velocity  along  this  boundary  will  not  be  equal  to  the  velocity  at  upstream 
infinity  since  some  inviscid  deceleration  will  have  occurred.  However,  as  long 
as  the  boundary  is  upstream  of  the  region  of  any  significant  viscous  or  shock 
phenomena,  the  total  pressure  on  this  boundary  will  be  equal  to  the  total  pres¬ 
sure  at  upstream  infinity.  Hence,  total  pressure  is  an  appropriate  boundary' 
condition.  In  addition  to  specifying  upstream  total  pressure,  it  is  necessary 
to  specify  the  inlet  flow  angle.  In  the  present  calculation  a  value  was  simply 
assumed.  A  more  accurate  specif icat ion  of  the  inlet  flow  angle  could  be  the 
embedding  of  a  viscous  region  surrounding  the  cascade  within  a  large  inviscid 
flow  field.  Zone  embedding  processes  have  been  discussed  in  some  detail  by 
McDonald  and  Briley  in  Ref.  35.  Obviously,  the  current  specif icat ion  is  a 
simple  one.  An  alternate  and  quite  feasible  specification  of  the  inlet  flow 
angle  could  be  based  upon  an  inviscid  solution  for  the  cascade  and  this  could 
be  improved  by  using  an  interaction  between  the  viscous  and  inviscid  analyses 
to  allow  viscous  phenomena  occurring  near  the  blades  to  affect  the  flow  angles 
at  the  upstream  boundary.  This  refinement  is  not  included  at  the  present  stage 
of  development  of  the  analysis.  The  third  condition  set  on  the  upstream  boundary 
concerns  the  density  and  several  possibilities  exist.  Either  the  first  or 
second  normal  derivatives  can  be  set  to  the  inviscid  value  or  be  set  to  zero. 

In  the  present  application  the  first  derivative  of  density  was  set  to  zero. 

The  downstream  boundary  was  treated  by  setting  a  constant  static  pressure 
and  by  setting  second  derivatives  of  both  velocity  components  equal  to  zero  at 
this  location.  In  the  present  application  the  static  pressure  was  set  as  the 
pressure  at  downstream  infinity,  and  hence  it  is  assumed  that  the  downstream 


35.  McDonald,  H.  and  W.  R.  Briley:  Computational  Fluid  Dynamic  Aspects  of 
Internal  Flows.  AIAA  Paper  No.  79-1445,  1979. 
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boundary  is  located  in  a  region  where  pressure  is  uniform.  Hnvever,  us  -.  it;; 
the  inlet  flow  angle  it  would  be  possible  to  calculate  the  pressure  along 
this  boundary  with  an  inviscid  analysis  and,  if  necessary,  to  interact  the 
viscous  and  inviscid  analyses  as  the  flow  progresses. 

Both  the  upstream  and  downstream  boundaries  have  boundary  conditions 
associated  with  them  which  are  nonlinear  functions  of  the  independent  variables. 
These  are  the  specifications  of  total  pressure  on  the  upstream  boundary’  and 
static  pressure  on  the  downstream  boundary.  These  nonlinear  boundary  condi¬ 
tions  are  linearized  in  the  same  manner  as  the  governing  equations,  and  then 
solved  implicitly  along  with  the  interior  point  equations. 

The  final  boundary  conditions  to  be  considered  are  the  conditions  along 
the  blade  surface.  Here  no-slip  and  no-through  flow  conditions  were  applied 
leading  to  a  specification  of  zero  velocity  on  the  surface.  An  inviscid  trans¬ 
verse  momentum  equation  was  applied  on  the  surface  leading  to  a  boundary’  con¬ 
dition  approximtion  of  2ero  transverse  pressure  gradient  being  applied. 
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IV. 


As  a  result  of  the  advances  made  in  the  coarc i rate  ‘■vsttr  fer.t rat; : on  - 
cess  it  becomes  practical  to  apply  the  Navier-S takes  calculation  procedure  to 
realistic  cascade  geometries.  The  geometry  chosen  is  the  so-called  Jose  Sana 
controlled  diffusion  cascade,  with  solidity,  c  =  c/ri,  equal  to  approximat e ly 
0.94,  an  inlet  flow  angle  of  35.  S°  and  an  exit  geometry  angle  of  0°.  For  the 
purpose  of  the  calculation,  the  trailing  edge  was  brought  to  a  cusp  to  allow  use 
of  a  ’C'-tvpe  coordinate  grid.  The  solidity  ratio  is  defined  as  the  ratio  of 
chord,  c,  to  distance  between  periodic  lines,  H.  Tne  flow’  conditions  chosen 
were  Reynolds  number  based  upon  chord  and  downstream  conditions  of  1.2  x  10°, 
which  gives  largely  turbulent  flow  over  the  cascade,  and  a  nominal  approach  Mach 
number  of  0.75.  Since  the  boundary  conditions  set  total  pressure  at  the  upstream 
boundary  and  static  pressure  at  the  downstream  boundary,  the  Mach  number  on  the 
upstream  boundary  is  not  set  a  priori  but  evolves  from  the  solution  and  is  con¬ 
trolled  by  the  ratio  of  these  quantities. 

The  calculations  were  carried  out  on  a  highly  stretched  grid  having  113 
pseudo-azimuthal  points  and  30  pseudo-radial  points;  i.e.,  referring  to  Fig.  10, 
113  pseudo-radial  lines  such  as  A  A1  and  30  pseudo-azimuthal  loops  A'E'D*  were 
used.  Grid  points  were  packed  to  obtain  high  resolution  both  in  the  vicinity  of 
the  blade  surface  and  in  the  vicinity  of  the  blade  leading  edge.  The  first  point 
off  the  surface  was  placed  at  0.00004  chords  from  the  surface;  in  contrast  the 
distance  between  grid  points  in  the  center  region  of  the  passage  reached  a  value 
of  0.06  chords.  Similarly  the  distance  between  grid  points  on  the  blade  surface 
in  the  vicinity  of  the  blade  stagnation  point  was  also  kept  small  having  a  value 
of  .0046  chords  at  the  blade  leading  edge.  In  contrast,  the  grid  spacing  at  the 
downstream  boundary  was  0.16  chords,  the  dowTistream  boundary  being  taken  two 
chords  downstream  of  the  blade  trailing  edge. 

The  calculations  were  initiated  by  assuming  uniform  flow  with  a  boundary 
layer  correction  on  the  airfoil  surface.  The  calculation  wras  run  first  as  a 
variable  viscosity  laminar  flow  and  after  the  basic  flow  pattern  began  to 
develop,  a  mixing  length  turbulence  model  was  initiated.  The  calculations  were 
run  using  a  time-scaling  technique  developed  by  McDonald  and  Briley  (Ref.  35), 
and  previously  used  for  cascade  flows  by  Shamroth,  Gibeling  and  McDonald 
(Ref.  13).  This  technique  in  conjunction  with  the  existing  LBI  Navier-Stokes 
calculation  procedure  (Ref.  14)  leads  to  a  very  efficient  numerical  procedure 
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which  c ;r.vr rg^d  rapidly,  Convergence  was  err airsed  :r:r.  the 
approximately  153  t  ime  steps  with  convergence  being  du  terr.ir.t_c  bv  r: 
the  tine  history  of  the  flow  at  selected  points  in  the  field  as  well 
tine  history  of  the  surface  pressure  distribution.  No  attempt  was  made  to 
optimize  convergence  so  it  may  be  possible  to  converge  the  solution  in  Itrss 
than  150  time  steps.  The  calculation  for  this  grid  required  apprerimat elv  15 
seconds  of  CDC  7600  time  per  time  step.  However,  the  code  used  to  implement 
the  algorithm  is  still  a  research  code  with  a  large  amount  of  overhead  computa¬ 
tion.  With  some  recoding  it  is  believed  run  time  per  tine  step  could  be 
reduced  substantially. 

Three  cases  were  calculated.  The  first  case  corresponds  closely  to  the 
inviscid  calculation  presented  in  Ref.  6  with  the  upstream  flow  angle  being 
30.8°  and  the  converged  upstream  Mach  number  being  0.72.  The  case  was  calculated 
using  second-order  artificial  dissipation  in  both  coordinate  directions  to  eli¬ 
minate  spatial  oscillations.  In  previous  subsonic  calculations  for  both  isola¬ 
ted  airfoils  (Refs.  32  and  33)  and  cascades  (Ref.  13)  reasonable  results  had 
been  obtained  with  the  dissipation  parameter  equal  to  0.5  in  both  the  pseudo- 
radial  and  pseudo-azimuthal  directions.  Although  these  values  appeared  reason¬ 
able  for  shock  free  flow,  the  one-dimensional  study  previously  described  indi¬ 
cates  a  lower  value  should  be  used  in  the  direction  normal  to  the  shock.  There¬ 
fore,  in  this  first  case  the  pseudo-azimuthal  dissipation  parameter  was  taken 
to  be  .05  with  an  intention  of  investigating  lower  values  (o  -  .025)  subsequently 
Since  the  model  Rayleigh  problem  focuses  upon  the  shock  capturing  aspects  of  the 
transonic  cascade  analysis,  it  does  not  give  guidance  as  to  a  choice  for  the 
dissipation  parameter  in  the  pseudo-radial  direction  where  normal  diffusion 
rather  than  axial  diffusion  is  a  major  contributor  to  the  momentum  balance. 
Therefore,  the  diffusion  parameter  ir.  this  direction  was  kept  at  the  previously 
used  value  of  0.5.  It  should  be  noted  that  with  the  small  pseudo-radial  grid 
spacing  in  the  vicinity  of  the  blade  boundary  (An/c  *  4  x  10  and  with  the 
small  transverse  velocities  in  this  region  (v/u  <  .05)  little  pseudo-radial 
artificial  dissipation  was  added  in  the  boundary  layer  region;  the  major  effect 
of  pseudo-radial  artificial  dissipation  is  in  the  cascade  center  region.  Al¬ 
though  this  appears  reasonable  particularly  since  the  shock  is  nearly  parallel 
to  the  pseudo-radial  direction,  the  role  of  this  artificial  dissipation  mecha¬ 
nism  must  still  be  investigated  further.  Such  a  study  is  envisaged  in  the 
near  future. 
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the  calculation  of  the  r  =  t  .  ..  . 
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tne  Vicinity  of  the  shock,  it  was  not  adecuate  in  the  leacir.c  c-d-e  ve-.-c^ 
mis  region  contained  very  large  spatial  oscillations  when  the  run  was  r.ade 
uncer  the  a  =  .05  conditions.  Therefore,  o  was  made  a  function  of  distance 
fret1,  the  leading  edge  being  taken  as  0.5  at  the  leading  edge  and  decreasing 
smoot-hly  to  a  value  of  0.05  by  trie  ten  per  cent  chord  station.  It  is  ar.tici- 
f “ v chat  the  nec. ess  it y  to  increase  z  in  the  leading  ecce  is  related  to  a 
lack  of  spatial  resolution  and  with  better  spatial  resolution  the  pseudo- 
<sZ^.utnal  dissipation  parameter  could  be  set  to  0.05  throughout  the  flow. 

Although  some  differences  between  viscous  and  inviscid  solutions  should 
exist,  they  should  be  consistent  with  the  inclusion  of  viscous  phenomena.  A 
comparison  between  the  inviscid  calculation  and  the  current  viscous  calculation 
witn  Re  ^  1.2  x  10  and  M  =  0.72  is  presented  in  Fig.  11.  .As  can  be  seen  in 
Fig.  11,  the  results  of  the  two  procedures  compare  quite  well  in  the  leading 
edge  region.  The  maximum  suction  peak  value  in  the  viscous  solution  is  slightlv 
higher  than  that  obtained  in  the  inviscid  solution;  however,  both  peaks  occur 
at  the  same  location.  Differences  in  the  trailing  edge  surface  pressure  distri¬ 
butions  are  consistent  with  the  inclusion  of  boundary  layer  development  in  the 
viscous  solution.  Although  a  full  quantitative  assessment  must  await  an  exten¬ 
sive  data  comparison  program,  the  results  shown  in  Fig.  11  demonstrate  the 
Kavier- Stokes  results  to  be  consistent  with  inviscid  results  at  high  chord 
Reynolds  numbers  and  low  incidence,  as  expected. 

The  second  and  third  cases  considered  serve  to  show  details  of  the  flow 
field  and  the  effect  of  the  dissipation  parameter,  o.  In  Case  II,  the  dissipa¬ 
tion  parameter  was  taken  to  be  0.5  in  both  the  pseudo- radial  and  pseudo-azimuthal 
directions.  The  ratio  of  upstream  total  pressure  to  downstream  static  pressure 
was  taken  as  1.48  and  the  upstream  Mach  number  of  the  converged  solution  was  ap¬ 
proximately  0.78.  The  third  case  considered  is  similar  to  Case  II  with  the 
dissipation  parameter  in  the  pseudo-radial  direction  being  taken  as  0.05. 

Rather  than  keep  the  ratio  of  upstream  total  to  downstream  static  pressure  con¬ 
sistent  with  that  of  Case  II,  the  downstream  static  pressure  was  varied  so  as 
to  match  the  inlet  Mach  number  of  0.78  at  the  upstream  boundary  in  both  cases. 

In  other  words  the  cases  were  run  at  the  same  upstream  total  conditions  with 
the  downstream  pressure  being  varied  to  obtain  the  same  upstream  Mach  number. 

A  Mach  number  of  0.78  was  obtained  in  this  case  when  the  ratio  of  upstream 


total  pressure  to  downstream  static  pressure  was  1.-0.  Th^s  is  ct  r.r  i  e  t  *-r.  t  with 
a  lover  loss  level  at  the  lower  value  c:  the  dissiy.-.t  icn  r^r^r.c  ter . 

Prediction  of  the  surface  pressure  distribution  for  the  two  values  of 
dissipation  parameter  is  shown  in  Fig.  12  where  the  effect  of  artificial  dis¬ 
sipation  is  evident.  Although  the  surface  pressure  distribution  rerr.2ir.ed 
smooth  when  the  dissipation  parameter  was  0.05,  spatial  oscillations  occurred 
in  the  densitv  field  away  from  the  wall  in  the  vicinity  of  the  shock  and  had  an 
amplitude  of  approximately  four  per  cent  of  the  mean  density  value.  These  os¬ 
cillations  were  not  present  in  the  model  problem  even  when  the  dissipation  para¬ 
meter  was  reduced  to  .025  and  their  appearance  here  is  indicative  of  the  more 
complex  flow  field  present  in  the  cascade.  Since  some  oscillations  appeared 
when  the  dissipation  parameter  was  taken  as  .05,  no  attempt  was  made  to  reduce 
this  value  to  .025  as  it  was  felt  that  this  would  increase  the  amplitude  of  the 
oscillation.  As  can  be  seen  in  the  figure,  setting  the  dissipation  parameter 
to  0.5  for  the  pseudo-azimuthal  (streanvise)  direction  gives  considerable 
non-physical  smoothing  in  the  region  of  the  shock  and  reduces  the  suction  peak 
from  a  value  of  -1.5  at  a  dissipation  parameter  of  .05  to  a  value  of  -1.16  at 
a  dissipation  parameter  of  .5.  Obviously,  the  use  of  the  dissipation  parameter 
equal  to  .5  for  the  azimuthal  direction  has  had  a  very  adverse  effect  upon  ac¬ 
curacy  for  this  calculation.  On  the  other  hand  the  0.05  dissipation  parameter 
gives  a  sharp  representation  of  the  pressure  rise.  The  s trearcvise  grid  spacing 
in  the  vicinity  of  the  shock  is  Ax/ c  .03  and  although  the  pressure  rise  is 
smeared  in  the  wall  region  due  to  the  shock  boundary  layer  interaction,  it  is 
very  sharp  in  regions  removed  from  the  boundary  layer.  For  example,  at  a 
distance  of  0.05c  above  the  suction  surface  the  pressure  coefficient  rises 
from  -1.44  to  -0.31  in  four  grid  points. 

Contours  of  Mach  number  for  the  two  cases  are  presented  in  Figs.  13  and  14. 
The  outer  boundaries  are  periodic  lines  and  the  plotting  routine  draws  the  co¬ 
ordinate  system  branch  cut  emanating  from  the  airfoil  trailing  edge.  For  the 
.5  dissipation  parameter  calculation.  Fig.  13,  the  maximum  Mach  number  reached 
is  approximately  1.17  and  the  Mach  number  decrease  cs  the  flow  proceeds  over 
the  suction  surface  is  very  gradual.  In  contrast  the  maximum  Mach  number  for 
the  .05  dissipation  parameter  calculation  is  1.56  and  a  rapid  deceleration 
resembling  a  shock  wave  occurs  downstream  of  this  maximum.  A  comparison  of 
Mach  number  contours  for  the  flow  field  above  the  suction  surface  shows  the 
strong  effect  of  dissipation  parameter  and  indicates  setting  the  dissipation 
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parameter  to  .5  is  inappropriate  with  this  n-eh  distribution  ,u.c  ds:.rif  ::r 
flows  having  shock  waves.  As  previously  discussed,  the  azimuthal  c  =  .05 
calculation  snowed  some  slight  spatial  oscillations  in  density  (  )  which 

may  require  further  attention,  however,  as  seen  in  the  figures,  realistic 
solutions  were  obtained.  It  should  be  noted  that  the  effect  of  artificial 
dissipation  was  considerably  less  on  the  pressure  surface  flow  field  and  in 
the  blade  wake  region  where  no  shock  occurred.  Velocity  contours  which  give 
the  magnitude  of  the  velocity  vector  are  presented  in  Fics.  15  and  16.  The 
results  here  are  consistent  with  the  Mach  number  results. 

Velocity  vector  plots  for  both  cases  are  given  in  Figs.  17  and  18.  For 
clarity,  only  the  velocity  vectors  at  selected  points  are  plotted.  The 
velocity  direction  is  specified  on  the  inflow'  boundary,  however,  the  flow 
angle  is  not  specified  on  the  outflow  boundary  and  the  observed  flow  turning 
results  from  the  calculation.  The  predicted  flow  angle  at  the  outflow  boundary 
(which  is  two  chords  downstream  of  the  trailing  edge  and  not  showTi  on  this 
figure)  is  approximately  6  degrees  in  both  cases.  The  predicted  turning  of  the 
flow  as  it  passes  through  the  cascade  is  evident  as  is  the  rapid  deceleration 
on  the  suction  surface.  As  shown  in  Figs.  17  and  18,  both  calculations  predict 
separation  to  occur  on  the  suction  surface.  Although  separation  is  possible, 
particularly  in  the  presence  of  a  shock,  this  result  requires  experimental 
confirmation.  It  is  possible  that  the  prediction  of  separation  may  be  influ¬ 
enced  by  the  simple  turbulence  model  used  and  a  comparison  with  data  may 
indicate  the  need  for  further  work  in  this  area.  The  vertical  arrows  occur¬ 
ring  at  the  airfoil  surface  are  data  points  having  zero  velocity  (no-slip 
conditions) . 

Contours  of  static  pressure  coefficient  are  presented  in  Figs.  19  and  20. 
Although  not  shown  on  this  plot,  the  stagnation  point  pressure  coefficients 
for  the  both  a  =  .5  and  a  =  . 05  cases  were  1.07.  This  is  somewhat  lower  than 
the  value  of  1.16  which  would  be  expected  from  inviscid  considerations  and  the 
discrepancv  may  be  due  to  viscous  effects,  numerical  truncation  or  a  failure  to 
resolve  the  Heimenz  layer.  The  differences  in  pressure  coefficient  between  the 
two  cases  are  again  more  evident  over  the  suction  surface.  The  ratio  of  outflow 
to  inflow  static  pressure  for  the  cases  is  1.03  and  1.1  respectively.  A  one¬ 
dimensional  inviscid  estimate  of  the  pressure  ratio  gives  pQUt/p^n  ~  1*28  and, 
therefore,  the  predicted  pressure  rise  is  less  than  the  inviscid  estimate.  This 
difference  represents  losses,  with  the  greater  losses  occurring  in  the  o  =  .5 
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caf-w,  the  case  having  larger  artificial  dissipation, 
pressure  coefficient,  (p  -  p  )  /  (1/2  c  u  ~)  ,  are  given  in  Fics.  21  an 
respectively.  These  low  total  pressure  regions  are  confined  to  the  b 
layer  and  the  wake. 


V. 


CONCLUDING  KL*  LARKS 

In  the  present  study,  several  questions  relevant  to  the  viscous 
cascade  problem  have  been^ posed  and  answered  by  considering  simple  flow 
situations.  Tnese  questions  focused*  upon  spatial  differencing,  specif icat ion 
of  boundary  conditions  and  use  of  artificial  dissipation  in  flows  containing 
shock  waves.  In  regard  to  the  first  of  these  problems,  the  node!  problem 
clearly  show*s  that  converged  solutions  can  be  obtained  using  spatial  central 
difference  representa t ions  in  both  subsonic  and  supersonic  portions  of  the 
flow.  A  study  of  boundary  conditions  in  a  simple  one-dimens icn2l  problem 
indicated  that  for  flow’s  having  subsonic  inflow  and  outflow  condi t ions,  specif i- 
cation  of  total  pressure  on  the  upstream  boundary  and  static  pressure  on  the 
downstream  boundary  is  satisfactory  and  physically  realistic.  Finally,  various 
methods  of  calculating  flow’s  w'ith  shock  waves  were  considered  by  solving  a 
one-dimens ional  flow  problem  with  heat  sources.  Among  the  methods  considered 
for  obtaining  stable  solutions  in  the  presence  of  shocks  were  second-order 
dissipation  methods,  fourth-order  dissipation  methods  and  pressure  damping 
methods.  The  results  obtained  in  this  model  problem  indicate  that  with  the 
envisaged  grids,  second-order  dissipation  with  a  relatively  low  dissipation 
parameter,  a  zz  .025,  is  a  suitable  method  for  use  in  transonic  flow  problems. 
Further,  it  was  found  preferable  to  add  dissipation  terms  containing  derivatives 
of  the  primitive  variables  p,  u,  w. 

Based  upon  these  preliminary  studies,  a  calculation  was  made  of  flow 
through  a  cascade  of  Jose  Sanz  airfoils.'  In  this  calculation  spatial  central 
differences  were  used  throughout  to  represent  spatial  derivatives;  boundary 
conditions  consisted  of  upstream  total  pressure,  total  temperature  and  flow 
angle  and  downstream  static  pressure  and  second-order  artificial  dissipation 
was  used  to  suppress  spatial  oscillations.  The  coordinate  system  used  was  a 
constructive  ’C’-type  system  containing  113  pseudo-azimuthal  points  and  30 
pseudo- radial  points.  The  grid  contained  high  resolution  in  the  vicinity  of 
the  blade  surface  and  the  blade  leading  edge  region.  A  calculation  made  for 
this  cascade  at  low  loading  predicted  a  surface  pressure  distribution  in  good 
agreement  with  inviscid  results;  the  differences  were  consistent  with  the  in¬ 
clusion  of  viscous  effects.  Two  additional  calculations  were  run  at  higher 
loading,  each  having  a  Reynolds  number  based  upon  chord  of  1.2  x  10^  and  a 
Mach  number  at  the  upstream  boundary  of  0.8.  In  these  ca 1 cul at  ions ,  the  dis- 
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sipation  parameter  in  the  pseudo-radial  direction  was  set  to  0.5.  In  the  first 
calculation  the  dissipation  parameter  in  the  pseudo-az imuthal  direction  was 
taken  to  be  0.5,  however,  in  the  second  calculation  this  value  was  taken  as 
0.05.  In  both  cases  an  embedded  supersonic  flow  region  appeared  over  the 
suction  surface.  In  the  case  with  dissipation  parameter  equal  to  .05  the 
supersonic  region  was  terminated  by  a  rapid  pressure  increase  representing  a 
shock  wave.  However,  in  the  0.5  dissipation  parameter  case  considerable 
smearing  of  the  results  occurred.  These  results  are  consistent  with  the  results 
obtained  for  the  model  one-dimensional  problem. 

The  results  obtained  in  the  0.05  dissipation  parameter  case  show  the 
expected  qualitative  features  of  transonic  cascade  flow  fields.  The  flow  ac¬ 
celeration  around  the  leading  edge  suction  suface,  the  flow  turning,  the  develop¬ 
ment  of  suet  ion  and  pressure  surface  boundary  layers  and  the  appearance  of  an 
embedded  shock  are  all  plainly  shovm  in  the  figures.  Of  particular  encourage¬ 
ment  is  the  attainment  of  results  for  this  complex  flow  field  using  a  dissipa¬ 
tion  parameter  criterion  of  .05  in  the  azimuthal  direction.  The  results  also 
show  that  setting  the  parameter  to  0.5  leads  to  considerable  shock  smearing. 

When  a  calculation  was  initiated  as  uniform  flow  with  a  boundary  layer 
correction,  convergence  was  obtained  within  150  time  steps;  no  specific  attempt 
has  been  made  to  optimize  convergence  and,  therefore,  it  is  likely  that  with 
improved  time  step  selection  the  number  of  time  steps  to  reach  convergence 
could  be  less  than  this  value.  The  calculations  presented  clearly  demonstrate 
the  practicality  of  efficiently  predicting  high  Reynolds  number  transonic 
cascade  flow  fields.  A  constructive  coordinate  system  has  been  modified  to 
allow  specification  of  realistic  cascades  and  calculations  have  converged 
relatively  quickly.  Although  further  efforts  are  called  for,  particularly  in 
the  area  of  quantitative  assessment  of  results,  turbulence  modelling  and 
further  investigation  of  artificial  dissipation,  the  results  obtained  to  date 
are  very  encouraging. 
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Fig.  10  -  Cascade  coordinate  grid  -  all  grid  points  not  included. 
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Fig.  14  -  Mach  number  contours,  a  c  ,05, 
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Fig.  20  -  Static  pressure  contours,  a  =  .05. 
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